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ABSTRACT 



This is a theoretical study to determine a practical design 
formula for the buckling load of simply supported rectangular plates 
under combined axial load and normal pressure. 

The theory and method used was based on Samuel Levy's 
integrated solution of the non-linear differential equations as de- 
rived by von Karman. The Icngth-to-width ratio (a/b) of the plates 
considered was extended up to a value of 4:1. For each lengtli-to- 
width ratio a set range of normal pressures was used in obt-'.ining 
the different buckling loads. 

The study confirmed the conclusion arrived at by Levy; that 
is, the normal pressure increases the buckling load. It was further 
observed and concluded that for low values of a/b (less than 1. 0) the 
buckling loads obtained in the manner described in Chapter III are 
impractical to attain. Yield failure will probably occur before the 
calculated critical load is reached. 



In general, however, the authors' recommended buckling load 
formula for a steel plate is of the form, 

T. .b 






K 



I2(l-iu ) 



where the value of K is scaled from Fig. XXXV. The basic limitation 
on the accuracy of this value of K is the approximated expression(s) of 
the deflection equation. The computed value of must be compared 
with the yield strength of the mate’ lal since the validity of the theory 
only holds true within the elastic region. 

Thesis Supervisor: J. Harvey Evans 

Title: Professor of Naval Architecture 
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I, INTRODUCTION 



Background 

The discussion of buckling of plates is usually based upon a 
linear differential equation derived under the assumption that the 
deflections of the plates are small In comparison with their thickness. 

In particular, present design practice on hull plating under combined 
normal pressure and axial load is often baaed on approximate critical 
buckling load formulas. The normal load has generally been relegated 
to minor roles. The reasons purported to supported this were very well 
founded, both from the standpoints of theory and practice. Samuel 
Levy and others [ 2] , (3) arrived at the conclusion that normal 
pressure increases the buckling load and thereby indicated that it 
would be conser vative design to neglect the effect of normal pressure. 
The same conclusion was stated by Bleich [ 1 ] . Bleich further made 
mention of the fact that deflection of ship’s hull plating usually does not 
exceed one-half the plate's thickness and therefore concluded that 
linearized theory could be applied. 

It is worthy to note at this point, however, that deflections 
exceeding half the plate's thickness are possible, especially in the 
case of simply supported plates. In those cases linear theory of 
plates no longer applies. As will be briefly explained In the sub- 
sequent paragraph, the problem becomes a non-linear stress problem 
since deflections of the order of magnitude of the plate thickness must 
be considered. 

Briefly stated, the nonlinearity of the differential equations has 
its origin in the fact that. In the case of large deflections, there is an 
interaction between the membrane stresses and the curvature of the 
plate. This interaction leads to non-linear terms in the eq' ations of 
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equilibrium of the plate elements. 

The complete differential equations of the problem were 
formulated by von Karman, who added those non-linear terms 
pertaining to the flexural rigidity of the plate. Timoshenko I 6] , 
Marguerre and Trefftz [ 8] derived the expressions for the strain 
energy of plates with large deflections. 

An attempt to develop the large deflection theory for rec- 
tangular plates under combined bending and longitudinal compression 
was made by Bengston [ 7] . His studies included plates with simply 
supported and clamped edges. However, Bengston Introduced, in 
the course of his analysis, certain arbitrary assumptions which in 
part contradict each othe r and it is very doubtful whether the results 
of his analysis can be considered as entirely correct. 

In a series of papers Samuel Levy gave exact theoretical 
solutions for rectangular plates with large deflections. The 
theories developed include simply supported and clamped plates. 

The problem of rectangular plates carrying longitudinal 
compression and normal pressure is of prime importance in the 
design of the hull plating of ships. It is closely related to the 
quest ion of buckling strength of plates. 

Statement of the Problem 

It is common knowledge that factors of safety are actually 
factors of ignorance which in some cases have remained unchanged 
for a period of time for no obvious reason. Factors of safety may 
be used to accoimt for one or more of the following; (a) material 
imperfections, (b) faulty construction practices, (c) material de- 
gradations due to corrosion, and (d) lack of ability to make rigorous 
calculations. The conservatlvehess of design as previously mentioned 
prompted the authors to ask in what way would conservative design 
affect the factor of safety; that is, does conservative design actually 
compound the factor of safety? If, for instance the safety factor is 
increased by only a few tenths of a percent, then one is justified in 
using linear theory. However, if the increase becomes considerable, 
the design may become uneconomical. 
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With the salient points stated in the previous paragraph in 
mind and also with the conclusion of Levy on the effect of normal 
pressure on the buckling load of plates, the authors came to the 
conclusion that there is a need for further investigations. First, 
these investigations afford the possibility of checking the accuracy 
of simpler approximate methods and formulas. Secondly, and 
actually the ultimate goal of this thesis, these Investigations would 
hopefully lead to a formulation of a simple design formula for plates 
subjected to the simultaneous action of normal and edge loads. 

A literature survey disclosed there is no simple formula which 
a designer can use with great facility. Fortunately, general solutions 
to the differential equations are available in the form of Fourier series. 
Though the analysis of these solutions is highly involved and the numeri- 
cal work for obtaining special solutions is very laborious, a method 
(Chapter III) has been devised which lends Itself to computer programming. 

This thesis is concerned mainly with plates simply supported 
on all edges and subjected to combined loadings of uniform normal 
pressure and axial load in one direction. It is further hoped that it 
will serve as a forerunner of later investigations based oh different 
boundary conditions. 
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n. THEORY 



Nomenclature (See Fig. I): 



a plate length in the x-direction. 

b plate width in the y-direction'. 

a = a/b 

h plate thickness, 

w plate deflection. 

X, y coordinate axes with origin at corner of plate. 

E Young's modulus. 

M Poisson's ratio. (Note: Tabulated results and figures are 

for steel with = 0. 300) 
p uniform normal pressure on plate, 

e average compressive strain at edges y = 0, b. 

P axial load on plate. 

F stress function. 

D = Eh^/l2 (l - /i^). Flexural rigidity of the plate. 
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Fig. I Schematic diagram of plate under axial load 
and normal pressure. 
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Fundamental Equations 

The fundamental partial differential equations governing the 
deformations of flat plates have been derived by von Karman. They 
are given by Samuel Levy [ 2] and Timoshenko | 6] . The exact math- 
ematical analysis to determine the buckling strength of a simply sup- 
ported flat plate under combined edge compression and normal loading 
Involves the integration of these equations; viz. , 
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The boimdary conditions for a simply supported rectangular 
plate to be satisfied by equations (l) and (2) are for deflection, w, 
and edge bending moments per unit length to be zero at the edges 
of the plate; viz. , 
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The complete solution is more fully explained in reference [ 2] 
and, without loss of generality, the deflection equation is here only 
reproduced in the following Fourier series. 



w 



E I 

m=l n* 1 



w^ „ sin 



mirx 



sln^p 



(5) 



where the undetermined constants, w__ must satisfy the relation(s) 

m, n '' 

expressed by equation (9) of the same reference. 

The average compressive stain, e, at the edges y = 0, b is 
computed from equation (11) of [2] as: 
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In this thesis, the deflection equation for large values of a/b 
(o » 1) Is approximated by the following expression. 



w ® w- - sin — sin ^ . sin sin ^ 

1,1 a b 3, 1 a b 



w,-. 



5, , * 



and for a/b - values In the vicinity of 1 and less thsin l, the approximated 
deflection equation is. 



+ w, , sin sin U + W, , Bln sin 
3, 1 a b 0,0 a b 



There are no special reasons for using two deflection equations 
.for different ranges of a values, except that equation (7) restricts the 
deflected shape of the plate to one sine wave across the width and a 
combination of four sine waves along its length. The errors involved 
by this approximation are expected to be less than five percent ( 2) , 

[ 3] . For this reason, equation (7) could be reduced to a fewer 
number of terms as a is reduced. This is evidenced by the deflection 
equation used for or = 3.0 in refeience { 2] . More specifically, the 
contribution of the higher-ordered deflection coefficients in the x- 
dlrcctlon becomes less significant as a is reduced from an n » 1. 

However, as a is further reduced the assumption of one sine 
wave across the width of the plate becomes incorrect since the 
higher-ordered deflection coefficients in the y- direction become 
increasingly significant and cannot be neglected unless the desired 
degree of accuracy is sacrificed. Thus, for such a values, equation 
(8) is used. As before, it is believed that the errors incurred by 
using a finite number of deflection coefficients v/ould give results 
within the desired accuracy (see Table 13 of [ 2] ). 
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The determination of the value of a applicable to equation (7) 
or (8) and the more general case of expanding the work of Levy et al 
on simply supported rectangular plates within the practical values of 
a/b are problems examined in this thesis. 
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m. PROCEDURE 



With some additional steps and explanatory comments incorporated, 
the procedure followed in solving for the deflection equations was essentially 
the same as that outlined in reference [ 3]. The steps were as follows: 

1. The family of four or six simultaneous cubic equations 

corresponding to the same number of unknown deflection coefficients 

was first obtained. This was done by solving for the coefficients b 

j p, q 

defined by equation (8) of reference [ 2] , and substituting them into 

equation (9) of the same reference. The results are shown in Appendix A. 

3 

2. Each of the resulting equations was divided by h . This step 
nondlmenslonalized the family of equations. 

3. V;»lues of w, ,/h, w- ,/h, w„ ,/h, etc. were estimated 
corresponding to chosen values of Pb/Eh^ and pb'*/Eh^ -- the non- 
dimensional axial load and normal pressure, respectively. This was 
the most delicate step since at higher values of Pb/Eh^ (normal pressure 
held constant) more than one solution was found possible (see { 3j also). 

In other words, the estimated w’s, when not properly chosen, could lead 
to solutions other than those corresponding to a continuous change in 
buckled form from zero axial load to the buckling load. 

It was therefore decided that a step-by-step method of estimating 
the w's be employed. In conjunction with the subsequent steps 4, 5, and 
6, this method proceeds as follows: 

a). Corresponding to three low values of Pb/Eh , including 

3 

Pb/Eh a 0, three sets of solutions were obtained using a reasonable 

3 

set of w-estimates. These values of Pb/Eh were hoped to lead 

to unique solutions. For these estimated w's, refe -ence [ 2] or 

[ 3] was used as a guide. 
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b). From (a) above, succeeding estimates of w's as 
axial load is increased, were based on the continued trend 

3 

of each w-vs-Pb/Eh curve; viz. , the slope and the rate of 
change of slope computed from the last three correct points 
were employed In an abridged Taylor series of the form, 

h h2 

w(xq + h) = w(xq) + w ' (xq) ^ + w* (xq) ^ , 

where, w(xq h) « next estimate of the particular w. 

w(xq) * last computed value of w. 

w * (xq) = elope obtained from the last two computed 
points. 



w' (xq) * rate of change of slope based on the last 
three computed values of w. 



4. The resulting cubic equations shown in Appendix A were 
linearized by expanding the right-hand side of each equation in a Taylor 
series In the neighborhood of the estimated values of ^/h, v/^ g/h, 
Wg ^/h, etc,, omitting higher ordered terms. 

5. Grout's method ( 5] was used in solving for the difference 
between the estimated w's and their improved values. 



6. Step 5 was repeated until such time tlat the calculated 
error sum was less than a test constant. The error sum was defined 
by the authors as the measure of closeness of successive approximations. 
In equation form, it is as follows; 



Error ** 



N 

y 

i« 1 



^i 



where, Xj represents the absolute value of the difference between an 
estimate w and its improved value. 

The test constant used for equation (7) was O.OOOl and for 
equation (8), O.OOOOOt. Why two different values of the test constant 



e 
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were used would become evident from a study of the results In 
Chapter IV. 

7. The ave age strain corresponding to each axial load was 
computed from equation (6). 

8. In the vicinity of the buckling load where the strain started 

to change more rapidly than the axial load, the w, . /h was made the 

3 ^ 

independent variable in place of Pb/Eh . The other deflection coef- 
ficients were computed in a similar fashion as before. 

4 4 

Values of the buckling load were obtained for values of pb /Eh 
equal to 2. 50, 7. 50, 12. 50, 18. 00, 24. 50, and 30. 00, and at values of 
a/b up to 4. 00. These were plotted versus a/b and analyzed for a 
possible comparison with Bryan’s classical solution. 

The whole procedure(s) had been programmed with the use of 
the IBM-7094. Details not bth 'rwise covered here may be found in 
the programs shown in Appendix B. 
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^ p 
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IV. RESULTS 



Tables I to XII are sample results of the computer programs. 
Tables I to VI were obtained using equation (8) for an a/b =1.00 and, 
the rest of the tables were for an a/b = 4.00 using equation (7). 

Figures II to XXXIII are the graphical results of the same 

programs. The curves drawn are the stress-strain curve analogy 

for plates with the non-dimensional axial load, Pb/Eh , as ordinate 

2 2 

and the non-dimensional average strain, eb /h , as abscissa. (Note: 
Not all results for the chosen a/b-values used were plotted. ) 

Tabulated values of the critical loads, rounded to the nearest 
hundredth, are shown in Tables XII and XIV. 

Figure XXXIV is a plot of the critical Pb/Eh^ versus a/b. The 
curves between zero a/b and the smallest value of a/b used have been 
extrapolated. The curves have been faired at points where the results 
of both equations (7) and (8) either become tangent or intersect. Re- 
finements of the curves in the vicinity of the cusps were made by ob- 
taining more results at the a/b-values concerned. These results are 
shown in Tables XIII and XIV. 
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Table I - Values of deflection coefficients for various values of axial compressive load in the 
x-direction, P, for simply supported rectangular plate, a/b = 1.00, p = 0.300. 
Normal pressure, p = 2. 50 Eh^/b"^. 
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Table II - Values of deflection coefficients for various values of axial compressive load in the 

X- direction, P, for simply supported rectangular plate, a/b = 1.00, u = 0.300. 

4 4 

Normal pressure, p=7.50Eh /b . 
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Table II (cont. ) - Values of deflection coefficients for various values of axial compressive load in 

the x-direction, P, for simply supported rectangular plate, a/b = 1.00, u » 0. 300. 

4 4 

Normal pressure, p =« 7. 50 Eh /b . 
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Values of deflection coefficients for various values of axial compressive load in 
the X- direction, P, for simply supported rectangulai’ plate, a/b= i. 00, /j“0. 300. 
Normal pressure, p*12. 50Eh /b . 
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Table IV - Values of deflection coefficients for various values of axial compressive load in the 

x-direction, P, for simply supported rectangular plate, a/b = 1.00, u = 0.300. 

4 4 

Normal pressure, p* 18. 00 Eh /b . 
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Table V- Values of deflection coefficients for various values of axial compressive load 

the X- direction, P, for simply supported rectangular plate, a/b*1.00, ^“0.300. 
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Normal pressure, p * 24. 50 Eh /b . 
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Table VII - Values of deflection coefficients for various values of axial 

compressive load In the x-direction, P, for simply supported 
rectangular plate, a/b » 4.00, u » 0.300. Normal pressure, 
p » 2. 50 Eh^lb^. 
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Table VIII - Values of deflection coefficients for various values of axial 

compressive load in the x-directlon, P, for simply supported 
rectangular plate, a/b = 4.00, /j » 0.300. Normal pressure, 
p = 7. 50 Eh^/b^. 
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Table IX - Values of deflection coefficients for various values of 

axial compressive load in the x- direction, P, for simply 
supported rectangular plate, a/b ■ 4. 00, p * 0. 300. 
Normal pressure, p *12.50 Eh /b . 
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T^bleX - 


Values of deflection coefficients for various values of axial 
compressive load in the x-dlrectlon, P, for simply supported 
rectangular plate, a/b * 4.00, u * 0.300. Normal pressure, 
p * 18.00 Eh^/b^. 
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Table XI - Values of deflection coefficients for various values of axial 

compressive load in the x-direction, P, for simply supported 
rectangular plate, a/b=4.00, u » 0.300. Normal pressure, 
p « 24. 50 Eh*^/b^. 
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Table XII - Values of deflection coefficients for various values of axial 



compressive load in the 
rectangular plate, a/b 



p * 30.00 Eh^/b'^. 



Pb 




^^3, 1 


e1? 


h 


h 


0. 


2. 286 


0. 643 


0. 500 


2. 318 


0. 662 


1.000 


2. 350 


0. 681 


1. 500 


2. 384 


0.701 


2.000 


2.418 


0. 720 


2. 500 


2.452 


0.739 


3.000 


2. 487 


0.759 


3. 500 


2. 522 


0.778 


4.000 


2. 558 


0. 797 


4. 500 


2. 594’* 


0. 816 


5.000 


2. 630 


0.835 



13. 331 


2.920 


0.996 


13. 337 


2.910 


0.993 


13. 339 


2.900 


0.990 


13. 339 


2.890 


0.987 


13. 337 


2.880 


0.984 



[-direction, P, for simply supported 
4. 00, u = 0. 300. Normal pressure. 
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Table XIII 



2 

Values of (Pb/Eh ) critical at various values of a/b 
4 4 

and pb /Eh . These are results using equation (8). 



a/b 






pb*^/Eh^ ** 






2.50 


7. 50 


12. 50 


18.00 


24. 50 


30.00 


0.60 


79. 51 


79. 60 


79.68 


79.78 


79.89 


79.98 


1.00 


27.56 


28.07 


28.59 


29.17 


29. 86 


30.46 


1.20 


16.89 


17.49 


18.06 


18.66 


19. 32 


19.86 


1.40 


11.33 


12.29 


13. 13 


13. 95 


14. 84 


15. 54 


1. 60 


7.90 


9.35 


10.48 


11.57 


12. 74 


13. 67 


1.70 


6.77 


8.41 


9.68 


10.93 


12.30 


13.43 


1.80 


5.93 


7.71 


9. 14 


10.60 


12.31 


13.86 


1.90 


5. 36 


7. 19 


8.82 


10.62 


13.20 


19. 35 


2.00 


4.96 


6.83 


8.72 


11.38 


17.03 


17. 65 


2. 10 


4.89 


6. 54 


8.96 


14.79 


15. 68 


16.38 
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Table XIV 



Values of a-t various values of a/b 

and pb^/Eh"^. These are results using equation (7). 



a/b 




* 


pb‘^/Eh'^ = 






2. 50 


7. 50 


12. 50 


18. 00 


24. 50 


30. 00 


4. 00 


4.20 


6. 82 


8. 68 


10. 26 


11. 98 


13. 34 


3. 95 


4. 16 


6. 78 


8. 77 


10. 33 


12. 03 


13. 37 


3. 90 


4. 16 


6. 72 


8. 87 


10. 41 


12.08 


13. 40 


3.85 


4.20 


6. 67 


8.97 


10. 50 


12. 14 


13.45 


3.80 


4. 30 


6. 61 


9. 08 


10. 58 


12. 21 


13. 50 


3. 60 


4. 12 


6.49 


9. 64 


11. 04 


12. 57 


13. 77 


3. 50 


4.03 


6. 49 


10.00 


11. 32 


12. 80 


13.96 


3.40 


3. 98 


6. 53 


10. 00 


11. 65 


13.06 


14. 19 


3. 20 


3.92 


6. 71 


9: 13 


12. 50 


13. 74 


14. 77 


3. 10 


3. 92 


6. 86 


8. 98 


13. 09 


14. 17 


15. 13 


3.08 


3.92 


6.89 


8. 97 


12. 17 


14. 26 


15.21 


3.04 


3. 92 


6. 97 


8. 96 


11. 69 


14.47 


15. 38 


3.00 


3. 91 


7. 05 


8. 96 


11. 42 


14. 70 


15. 57 


2. 95 


3.91 


7. 17 


8. 98 


11. 19 


15. 03 


15.82 


2. 90 


3. 92 


7. 30 


9.03 


11. 05 


13. 98 


16. 11 


2. 85 


3. 94 


7.44 


9. 10 


10. 97 


13.40 


16. 47 


2.80 


3. 99 


7. 60 


9. 19 


10. 93 


13.08 


15. 19 


2. 60 


4. 02 


8.44 


9. 79 


11. 17 


12. 71 


13. 98 


2. 57 


4.04 


6. 89 


9. 91 


11. 25 


12. 74 


13.95 


2. 50 


4. 10 


6. 50 


10. 24 


11. 50 


12. 87 


13. 97 


2. 40 


4.20 


6. 33 


10. 83 


11. 97 


13. 20 


14. 18 


2. 30 


4. 30 


6. 31 


11. 58 


12. 60 


13. 72 


14. 61 


2. 25 


4. 36 


6. 33 


9. 15 


12. 99 


14.06 


14. 90 


2. 20 


4. 48 


6. 39 


8. 81 


13. 44 


14. 45 


15. 26 


2. 15 


4. 53 


6. 47 


8. 64 


11. 64 


14. 91 


15. 67 


2. 10 


4. 69 


6. 57 


8. 56 


10.93 


15. 43 


16. 16 


2.05 


4. 82 


6. 70 


8. 56 


10. 61 


13. 12 




2.00 


4. 97 


6. 86 


8. 61 


10.45 


12. 55 


14. 30 


1. 90 




7. 27 


8. 69 


10. 45 


12. 15 


13.48 


1. 70 


6. 87 


8. 68 


10.06 


11.36 


12.71 


13. 74 
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V. DISCUSSION OF RESULTS 



The most remarkable part of the results was the way the 
critical buckling load behaved with respect to varying values of 
a/b (Fig. XXXIV). A< can be seen in the corrolary figure (page 64) 
to Fig. XXXIV, there is a striking resemblance to B yan's 
buckling load solution of rectangular plates solely under the 
action of edge loadings. It may be recalled that the wave-like 
pattern of B? yan's curve of K versus a/b was due to the nature of 
the K-formula and to the minimizing procedure adopted in obtaining 
the critical load for any given value of a/b. Bryan's K-formula is 
a function of the plate's dimensions and of the number of half- sine 
waves the plate takes after buckling. Thus, it would seem from the 
foregoing statements to presuppose that the patterns obtained in this 
thesis are also due in part to the number of half-sine waves. 

At first glance the results would give no indication as to the 
number of waves the plate takes at the incipience of buckling. But, 
checks had been made by the authors on the relationships between 
deflections at the midwidth of the plate and at various distances along 
its length. In other wor ds plots were made of the deflection ratio, 
w/h, versus the distance ratio, x/a, at y = b/2. These were done 
for different (a/b)'s at the same normal pressure. For example, it 
was found that for a p = 18.00 Eh~/b there were five buckles at 
a/b =« 2.40. (The two half-sine waves resulting from the initial 
general downward deflection of the plate due to normal pressure 
were not included in the count. ) When these checks were expanded to 
a wide range of values of a/b, it was fovmd that points lying on the 
same trough of a p- curve gave the same number of buckles. 

Since only representative tables are given (Tables I to XII) it 
was deemrd unnecessary to insert all these checks in this thesis. If 
verifications are to be made on this particular point the authors suggest 
going through the programs shown in Appendix B. 
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There Is, however, no positive explanation as to why only 
an odd number of buckles is possible. There are either one, three, 
five, etc., buckles. An even number of buckles, which are obtainable 
in the case of Bryan's, seems entirely out of the question. On this 
particular aspect, a possible explanation that the authors believed 
could have happened was in the manner the deflection equations were 
obtained. The approximated deflection equations v/ere follow-ups of 
those used in the references [ 2] , [ 3) , [ 4],. The equations only have 
the odd-numbered, subscripted deflection coefficients which, in the 
final analysis, resulted in an odd function for the w- equation. At this 
point, there is no way of predicting what the results will be if the even- 
numbered subscripts of the deflection coefficients are also taken into 
consideration. 

In reference to the reliability of the results, it is to be noted 
that the basic theory is vali 1 only within the range of Hooke's law. 

No allov/ance was made in the original differential equations on the 
change of Young's modulus when the axial load exceeds the propor- 
tional limit. Hence, all the results of this thesis presupposed that 
the material had been within the elastic re^on at all times. Experi- 
mental results are not available on the case of plates with the same 
loading and boundary conditions such that, definite Identifications 
cannot be made as to whether the results are within the elastic 
region or not. 

When a further study is made of Fig. XXXIV for a/b less than 
or equal to one, there is a very wide discrepancy between the calcu- 
lated critical loads and Bryan's that it is quite doubtful whether the 
expected critical loads are below the yield strength of steel for a 
majority of plate characteristics. 

The following computed results are hoped to illustrate the point: 

Consider a plate with a b/h *= 104. 7 and acted on by a pressure, 
p = 12. 50 Eh^/b"^ = 3. 11 psi. The b/h may correspond to a 14. 02- 
pound plate with b = 36 inches, and the p may correspond to a hydro- 
static head of 7 feet. 
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Case 1 



Case 2 



Case 3 



a/b * 1, 50 



a/b = 1.00 



a/b = 0. 60 



Ks 12.8 



K s 28 



K a 78 



a a 31800 psi 
c 



a « 69000 pci a » 192000 pai 

c* c 



If mild steel or high-tensile steel (HTS) were used with yield 
stresses of 33000 and 47000 psi, respectively, it is apparent that 
cases 2 and 3 exceed their corresponding yield strengths. Hence, 
the theory fails for these cases. If HY80 or HYIOO were used then 
theory is found to fail only for case 3. At any rate, the fact remains 
that the basic theory does not differentiate between the clastic and 
inelastic regions. A consequent conclusion, therefore, is the need 
for a check oh the yield strength whenever the formulated of this 
thesis is used. 

Except for the restrictions mentioned in the pre’/ious paragraph, 
it is obvious that the results obtained confirmed the conclusion of Levy 
on the effect of normal pressure on the buckling load. But, it is also 
evident that to neglect the normal pressure in predicting buckling strength 
is unwise design from an economic point of view. 

Lastly, the computer programs are limited in scope, although 
they worked propei ly in providing all the needed results of this tliesis. 
They were the consequent developments of the approximated deflections 
equations. A more desirable program would be to start with the 
equations shown in reference ( 2] . 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



ConcluBions 

1. Normal pressure always increases the buckling strength 
of rectangular flat plates. However, the buckling loads obtained by 
means of Levy's solution are not conclusive values; that is, the 
buckling load may or may not be within the proportional limit. In 
other words, theory does not differentiate between the elastic and 
inelastic regions. 

2. Plotted values of critical loads against the ratio, a/b, 
for a constant value of normal pressure exhibits the same trend 
and characteristics as the plot of Bryan's formula for zero 
pressure. 

3. To neglect the effect of normal pressure in the determination 
of buckling load is unwise design from an economic standpoint. 

4. In general, the buckling load formula can be written in the 

form, 

= .7-1 K 

where, this time, K is not only .9 function of plate geometry but also 
a function of the normal pressure and Poisson's ratio as well. 

P ecommendations 

1. A better representative plot of K versus a/b is needed by 
including the even-numbered subscripted deflection coefficients of 
Levy's solution. 

2. In consonance with the first recommendation, a ^ evised 
computer prog- am is desirable, possibly to start right at the series 
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forms of the equations shown in reference [2] . This step would 
eliminate using approximate deflection equations. 

3. It is recommended that a study be made of the most 
feasible experimental set-up to simulate the boundary conditions 
stated in this thesis; and hence, the ensuing experimental 
verifications. 
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APPENDIX A 

DEFLECTION EQUATIONS 



The six (6) simultaneous cubic equations involved in the solution 
of the deflection equation. 



w = j sin^sin-^+ g sin^sin-j^ 






^ 1 sin — T — sin Wo o sin 
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+ w, ^ sini^sin-^ +W,, , gin ^ sin ^ 
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°1 pb_^ + 4 h^ 
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corresponding significance of each subscripted a is shown at the end of 
this Appendix. 
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* 3*3, 1 *3. 3 <“*1 * 25^2 * 8‘»4 * *®®5> '“ 1 . 1 *5, 1 

- 2(25a2 ■•' 441a^ + 256ag)Wj 3^5 1^ ’?29(a^ 

■♦■ an) w, , w„ nW, (. + (256a« + 2401a. •♦• 4096a- +2401a,n)w? «w 



‘8' "1, 1 "3, 3" 1. 5 



12'”1. 3"5, 1 



+ (4a, + 225H2 +256.3 + aj)W 3 _jW 5 _, + 3 *3, 1 "l, 5 

- 9(64.3+ 441.,+ 169a,3)w,_3W3_3W,_ 3 - 6(64.3 + 441.3 

+ 169.,2) * 3 _ , * 3 _ 3 w,_ 5 - 18(25.3 + 64.3 + 9 . 3 ) , v» 3 _ 3 , 

+ 1296 . 3 w|_ 3W,_ 5+ 81(16.3 + .,+ 81.3+81.3) g Wg_ , 

- (25a, + 576.3 + 8281a, 3>W3 ,w^ 3 - 2 ( 25.3 ^^ 6.3 

+ 8281a, «)w, _w, ,w, , + (1552a_+ 28561.,,, + 28561a, «)w, w_ , 

1, O 1, D D, 1 O 16 1, D t>, 1 

+ (aj^ + 625a2) Wg j. 

The four (4) simultaneous equations involved in the solution of 
the deflection equation, 

V/ a w. . sin sin n~- + Wn - sin sin ^ 

1, 1 a D 3, 1 a D 



+ w- , sin 
0 . 1 



5ttx 



sin 



15 " 



+ w,_, sin 



Ittx 

a 



sin . 



are as follows: 



pfi). * r “ g g g j I ^ '*' c*^'^ ^ j*s)a - 

I, s <.*^^“‘ •e“>s ‘,.c) f,.£»,.e ' 

rr - j60-:s * J*t» * J.«')C 4 . J jvj jWjXIW - 

I .»* .i" ^e*"' ‘ t*'-' * " e .r- i .e*8 .i’ * 

. ''■'p f*^'r . %-(_ fii'CE * . *■ nStDC. - 

j ,> c 1 ^ I p ♦ 

4 f. 38 C) + , _j« j _jW , ,- (g* ■• 

8 .j" 1 .8* t. .i'"e'*'''’ * ; .o'" I ,£■* 'd* * 8*®*®* s'*®® ‘ i'*> * 

gt;i.« 4 +f,«i‘* +e«»«)8 - 

I .e** 8 .8* 1 .8* *8“® ■’ c"® * ' - 8 , 1"6 .e*"! .e*<st*®'-'' ^ 

J ,W p ®W (gBlS 4 ^ati 4 ^*4 ^ B^I)I8 4 ^ jWg g-WgBCeSl 4 

„pS\*. 4 .»6S)S - - ?• , ,»•(. ,.-l»r.8 4 tsra * ,e6S) - 



j. jS> 



lo noi*wlot •di ffl l:9vld\ni iv'r-jtsaps a ucs ns it) ri 0 t'i »dT 
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0 £ 1^0 » M 1«>4 






' d 



^ t’ <7 *» X t 

— nxB , -V/ -ft- f.'.?* - nx8 . •* ♦ 
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0 



- -V ^ *1, 1 • ^ Ir *1. 1 * “i, 1 

- 3.jwf_jW3_,+ (-Ia, + 9a2+ IBaj + a^) w,_ j ^ 

+ (3a, + 64a3>w^_ j «5, 1 ' * ^®*4> ”'1 1 *7, 1 * <‘‘=1 

+ 25ag+ 81a^+16ag)Wj 1^5 j + (4aj ■•• 49ag 256ag 
4 81ag)w,_,w2_, - 6(a,4 3a^)w,_,W3_,W3 , - 2(3a, 

+ 64ag)Wj jWg jW^ j + ( 6 a^+ ^^^^3 225a,^ 

4 9ag) W3,,W5 _,w,^,. 



256^1 pb^ 



7 ^ ^ w . ii4 Pb - a 

Stt ^ 3(l-,u^) ®5 ^ iT^ 3, 1 1 1, 1 

+ (4aj ■>■ Oag + I6ag+ a^)Wj^ jWg ^ - 3(aj + 3a^) w“ , w. 

+ (aj + Blag) Wg j + (4aj + 225ag + 256ag + ag) 

2 2 
+ 3(a^ + 192ag) Wg j + (4a^ + 441ag + 625a^+ 16a^)Wg j 

4 2(3a,4 64a3)w,_,W3_,W3_, - 2(3a, 426a^) w,_ , Wj_ , w,_ , 

+ 3(2a, +48a„ + 75aa + 3a«) w, , w_ 

1 0 4 o 1, 1 5, 1 /, 1 



pb'^ _ 

Sd-Ai*') ”7 



2 : i:; ^^5. 1 2 m ^ 5,1 



400 Pb 



- 3(aj + 3a^)w^ 1 ^3 l ^ ISa^) j Wg ^ 

2 2 

- .(3aj + 64ag)Wj ^ w.^ ^ <3aj + 64ag) ^Wg ^ + (4aj + 225ag 
+ 256ag+ 8g) Wg jWg j + (a^+ 62582 ) Wg 1 ^ (4aj * 122582 
+ 1.298ag ■♦• a^)Wg j ^ + S(2aj + 48ag + 75a,j + 3ag) ^ Wg ^ ^ 

4 6(a,4 192a3)w3_,W3_,w,_,. 






8 , . . ^ d'i 3J 



w ,.J— ^ i! 






•*' 
• yl> 



ru-i)g 



n? 



f ,J’" -I*^' *-'• " ['•* ■* I .c* I ,t* r ^ ■ 

jSJ» + , ,.« J jw<,s3S * ,,j) - - 






\st + ^p '8 ^ ^ J I" 



‘^ ► , w{, tf 4 4 



^*CSS *- gfi^tr + ^Bd) 4 ^ ^ ^ffi^shd 4 



■i,-.'''i,s'''i,t''<o'“ " 



I .1 



V. cN »4f ^ . dq I*” '■" 

^ • " , u'W' . 7 sr T< — ~ * c.^ IfT “5 * “ 



^ jf 1 !~ ■ f ,< 
It 



f 4 



»/ 



' .l’" 't*® * J*'® " t ,* ( .!**•* * *.S* *l‘' 



! j 1 .e"*t‘ *£*'■’* * i**'* ‘t“ >+ t ’ !-'' *■ 

g 

W . ^ W ( “* • ^S£ 4 ^,8/i* > 4 »| ) 4 _VI ,W f^cSC \ ffijS 4 

f m** , C. A . ,» (.e*'S 4 - ,p»iC ' , *.w , , 4 ,rffc)S * 

. . kW i^et 4 fiPV 4 ,aH ^ ,«U>J 4 

I ,' J ,G l ,j o o f 

^ d^ OQ. 



r ” r T — ~ r “'•• f 

[ * n I - ^ r , , ^e, 



*dq 

in ■ 

T,' 



- 0 



s s 

. •/ , ( m 8 a - ,SfV 4 4 ,H*|) •;- -W . , w (.»8 4 .ft)' 

i 4*. i , f *. 4 6 • E ,8 t J r I 
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256a 



0 w • 









1 pb’ ^ 4 h“ 784 Pb 

2 _ , . «._ » .2 



- ( 3 a^ + 64 aj^)Wj j ^5 ^ '*' '^S 1 ''^7 1 

- (3aj + 2f^a^)wJ ^ Wg ^ + (4 rj + 441ag ■♦• 625a^ + I6a^) Wg j ^ 

+ 3(aj + 192ag)Wg ^Wg ^ 4- (4aj + I225a2 + I296ag + a^^) w^ 1 ''*'^7 1 
+ (hj + 240182 ) 1 3^2aj + 48ag + 75a^ + Sag) ^ Wg 



=» O' 



The subscripted a's have the following significance: 

> 2 2^ 

Sg » a l(Qa^ + 1) 



Sg ® l/a‘ 



aio = a ^/(ISa^ + 1) 



8g « ce'^ f{a^ + 1) 
a^ » f(a^ + 4) 



ajj = a^H2'M^ + 1)' 



8j2 ” H^o/^ + 9) 



2 



2 2 * 
8g * a /(a + 9) 



8j3 ^ 4) 



2 // 2 ^ ,„.2 
Sg « O' /{a +16) 



aj^ + 36) 



a^ = + 25)^ 



8jg s? la^ + 49) 



ag = a^(4a^ + 1) 



where or » a/b. 












w 
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M3888-3^01»FMS»DEPUGU»1 ’500^0 BR I LLANTES-L I VADAS 

SOLUTION OF THE DEFLECTION EQUATION 

¥ = + V-3,i sv-n-M 

+V5 , sm 121 + ^7.1 ■— 

■’) cx b ^ b 

THIS METHOD IS SET UP FUR A CONTINUOUS 
CHANGE IN axial LOAD FROM ZERO TO A VALUE '.'ITHIN 
THE VICINITY OF ThE CRITICAL BUCKLING LOAQ. ZERO 
LATERAL LOAD IS EXCLUDED. 

MAIN PROGRAM 



XEQ 
LABEL 
L IST8 
FORTRAN 
READ 1» IND 

1 FORMAT (IS) 

PRINT 2 ' 

2 FORMAT (82H1 PNORM AXIAL W(l»l)/H W(3,l)/H W(5»l)/H W(7. 

ID/H strain error ) 

DO 44 JOB = 1 , IND 

READ 3. U, ALPHA* PNORM, LOOP* LAP 

3 FORMAT (3F8. 0*218) 

DIMENSION 31(4), B2(4)* B3(4) 

CALL CONST (U* ALPHA* PNORM* A1 * A2* A3* A4* AS* A6* A7* A8 * A9 * 

1 AlO* All* M2* A13* A1'4* A15* A16* A17* B1 * B2 * B3 ) 

INITIAL ESTIMATES OF .M I ) ' S 
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DIMENSION W(4) 

W ( 1 ) = 0.400 
W(2) = 0.300 
W(3) = 0.200 
W ( 4 ) = 0.100 

C the f^OLLOWlNG IS ThE ROUTINE FOR 

C CALCULATING THE DEFLECTION COEFFICIENTS FOR 

C CHOSEN values OF AXIAL LOAD AND NORMAL LOAD 

DIMENSION A(4,4), C(4)» X(4), Z(3), S(3»4) 

AX = 0. 

MG - 1 

DO LL = 1,3 

4 CALL SESAME 

CALL CROUT (A, C, X, M) 

GO TO ( 5 , 42 ) , M 

5 ERROR = 0. 

DO 6 J = 1 ,4 

6 ERROR = ERROR + ARSF(X(J)) 

IF ( ERROR - 0.0001 ) 9, 9,7 

7 DO 8 1=1,4 

8 W(I) = W(I) + X(I) 

GO TO 4 

9 CALL CPRNT (PNORM, AX, W, ERROR, Al) 

CALL STORE (AX, V.', Z, S, LL) 

AX = AX + 0.5 
10 CONTINUE 
DEL = -^.5 
KO = 1 



Ml - A* Q^lll* l ‘f*l I 

null 4 1 hS 1^1 M II I - J * ^ A. * • 

• 

••Of ^ -I . €• d 
»*•*» • «♦^> II iiif Cilll|l| | 




c 



>I 4 » M 



IP 



• «!«• 
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22 


GO TO 


( 23, 


25 


23 


del 


= 


0.05 






DO 24 


I = 


1 , 




W ( I ) 


= 


S(2 


» I ) 


24 


S ( 3 , 


I ) 


= W ( I ) 




AX = 


Z ( 2 ) 






Z(3 ) 


= 


AX 






AX = 


AX 


DEL 




KO = 


2 








KK = 


0 








CALL 


EXCH 


(Z, 




GO TO 


13 




25 


DEL 


= 


0.01 






DO 26 


I = 


1 » 




W( I ) 


= 


S(2 


» 1 ) 


26 


S ( 3 » 


I ) 


= W (I ) 




AX = 


Z( 2) 






Z ( 3 ) 


= 


AX 






AX = 


AX + 


DEL 




<0 = 


3 








KK = 


0 








CALL 


EXCH 


(Z, 




GO " 


0 


13 





THE following IS THE ROUTINE WHEN IN THE 
VICINITY OF THF BUCKLING LOAD. W(l) IS MADE 
THF independent VARIABLE IN PLACE OF THE AXIAL 
LOAD. 



21 DO 28 I 



1 ,A 







V. 

-i §5' 



I 

•» 
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31 FRROR = 

DO 32 J = 1,4 

32 ERROR = ERROR + APSF(X(J)) 

IF (ERROR - 0.0001) 36, 36, 33 

33 KK = KK + 1 

IF (KK - LOOP) 34, 34, 40 

34 DO 35 1=2,4 

35 W(I) = W(I) + X(I) 

AX = AX + X ( 1 ) 

GO TO 30 

36 CALL CPRMT (PNOPM, AX, W, ERROR, Al) 

IF (AX - Z(2 ) ) 44, 44, 37 

37 Z ( 3 ) = AX 

DO 38 I = 1,4 

38 S(3, I ) = ^( I ) 

DLl = 7(3) - Z(2) 

39 CONTINUE 
GO TO 44 

40 PRINT 41 

41 FORMAT (40H INFINITE LOOP OR NEEDS MORE ITERATION) 

CALL CPRNT (PNORfS AX, W, ERROR, Al) 

GO TO 44 

42 PRINT A3 

43 FORMAT (19H CROUT WENT CRAZY) 

44 CONTINUE 
CALL EXIT 

COM*^ON W, AX, B1 » B2, B3 , Al, A2, A^, A4, A3, A6, A7, AP* Ac. 
I All, AT 2, AT 3, A) 4, A13, A i 6 , \ 1 , M-A , A, C 



AID , 



MIQ > 
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END 



SUBROUT I NES 



C SUBROUTINE FOR SOLVING THE LINEARIZED SYST'^m 

C OF EQUATIONS BY THE GROUT REDOX METHOD 

» LABEL 

» LISTS 

* FORTRAM 



SUBROUTINE GROUT (A, G» X, M) 

DIMENSION A(4>/)» G(4)» AA(4,4)» GG(4)» X(4) 
DO 100 I = 1»4 

100 AA ( I , 1 ) = A ( I , 1 ) 

DO 101 J = 2»4 

101 AA ( 1 , J ) = A ( 1 » J ) /A ( 1 , 1 ) 

DO 102 I = 2>4 

DO 102 J = 2»4 

102 AA ( I » J ) = 0. 



J = 2 

103 II = J 

DO 105 I = I I ,4 
LIMl = J - 1 
DO 104 K = 1,LIM1 

104 AA(I»J) = AA(I,J) + (AA(I,K) * AA(K»J)) 

105 AA(I»J) = A(I,J) - AA(I,J) 

IF (AA( J,J) ) 106, 116, 106 

106 IF (4 - J) 110, 110, 107 

107 I = J 



m *■ ^ 




I 4MI 

'll i ^ * - 



Ml 4 till * • - 

- MTUl^ 04*.*® ^ 

... . •♦ «|«i Y^-**^* 



I 
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J = J + 


1 




JJ = J 






DO lOQ 


J = JJ,4 




LIM2 = 


I - 1 




DO 108 


K = 1,LI^^2 


lOfl 


AA ( I , J ) 


= AA(I,J) + (AA(I, 




AA ( I , J ) 


= A ( I , J) - AA( I , J) 


109 


AA { ! , J ) 


= AA( I ,J)/AA( 1,1) 




J = I + 


1 




GO TO 103 


110 


DO 111 


1 = 1,4 




X ( I ) = 


0. 


111 


CC ( I ) = 


0. 




CC ( 1 ) = 


C(1)/AA(1,1 ) 




DO 113 


I = 2,4 




LI M3 = 


I - 1 




DO 112 


K = 1 , L’M3 



11? CC(I) = CC(I) + (AA(I»K) * CC(0) 
CC ( I ) = C ( I ) - CC ( I ) 

113 CC(I ) = CC( I )/AA( 1,1) 

X ( 4 ) = CC ( 4 ) 

DO 115 1=1,3 

II = 4 - I 
LIM4 = II + 1 
DO 114 < = LIM4.4 

114 X(II) = X(II) + (AA(II,K) * X(K)) 

115 X(II) = CC(II) - X(in 

M = 1 




tw 



I • 

fi-il i •' • • 



• I • • 



fi - 1 1FI 



I «4 
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GO TO 118 

116 PRINT 117 

117 FORMAT (28H SINGULARITY - NO SOLUTION) 

M = 2 

118 RETURN 
FND 

C SUBROUTINE FOR CALCULATING THF NON-D I ME VS I C N A L 

C STRAIN AND THEN PRINTING ALL THE NEEDED RESULTS 

» LABEL 

» LIST8 

» FORTRAN 

subroutine cprnt (pnorm, ax, w, error, A1) 

DIMENSION W(4) 

SUM = 0. 

DO 200 I J = 1 ,4 
XY = 2 * IJ - 1 

200 SUM = SUM + (XY * W(IJ))**2 

STRAIN = AX + 3.1416**2 * SUM/(8. * Al) 

PRINT 201, PNORM, AX, (W(I), 1=1,4), STRAIN, ERROR 

201 FORMAT (F8.2, Fll.5, 4F10.3, Fll.5, F12.5) 

RETURN 

PND 

C SUBROUTINE FOR STORAGE OF THREE SETS OF VALUES 

C OF AXIAL LOAD /'ND THF FOUR DEFLECTION COEFFICIENTS. 

* LABEL 

» LISTS 

* fortoam 



SUBROUTINE STORE (AX, W, Z, S, LL) 



I I 



I 




lYl* 



I 






^ • « 



M mU 



« 



I 



4 



i§m§ 



I I 






|>r l.a«^ l — 



).l 



IJi« 



• ^11 i( * '■ 

,.,ii <•» • •> 

|M • ^ 11 ^ «•••!• 

, l!>»i<h. .-Mil— IT |<0 

, 1 . 1 •• ^ 







l«l 



V 
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DIMENSION W(4), Z(3), 5(3,4) 







Z ( LL ) = 


AX 






DO 203 


II 




203 


5 (LL , I ) 


= W ( I ) 






RETURN 








END 




c 






THIS subroutine 


c 






OF values from 






LABEL 




* 




LISTP 








FORTRAN 








SUBROUTINE EXCH (Z, S) 






D I MENSION Z ( 3 ) , S ( 3 ,4 ) 






DO 2^4 


1 = 1,4 






DO 204 


J = 1 , 2 






JA = J 


+ 1 




204 


S( J, I ) 


= 5( JA, I ) 






DO 205 


IJ = 1,2 






JP = IJ 


+ 1 




205 


Z ( I J ) = 


Z ( JB) 






DO 206 


J I = 1 




206 


S ( 3 , J I ) 


= 0. 






Z(3 ) = 


0. 






R'^TURN 





'NO 



subroutine for covputino all the necessary 
constants in TH'^ '-^AIN program. 



LAB-L 






f 






i :: i>\ - ^ 

- 

• 04 * 




-« V. 







'fm «» 



. I • 

< lltl 
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LIST8 

FORTRAN 

SUBROUTINE CONST (U* ALPHA » PNOR^S Al» A2» A3» A4 , A5» A6, A7, A8> 
1A9> A10» A11> A12» A13> A14» A15» \1A> A17» B1» ^3) 

DIMENSION Bl(4)» B2(4), B3(4)» P(9) 

P(l) = AL^HA**2 
P(2) = l./P( 1) 

DO 1 = 1,7 

Y1 = I**2 

11=1+2 

300 P(I1) = P(1)/(P(1) + Yl)**2 

A1 = P( 1 ) 

A2 = P ( 1 ) + P ( 2 ) 

A3 = P( 1 ) + 81. * P(2 ) 

A4 = D( 1 ) + 625 . * P(2 ) 



A5 


= 


P( 1 


) 


+ 24( 


1 . 


* P ( 2 ) 
















A6 


= 


? . 


* 


P( 1 ) 


JL 


64 . 


* P ( 3 ) 
















A7 


= 


3. 


* 


P( 1 ) 




576. 


* P(3) 
















A P 


= 


3 . 


* 


P ( 1 ) 


+ 


9. 


rr P ( 4 ) 
















A9 


= 


3. 


* 


P ( 1 ) 


+ 


25. 


* p ( 4 ) 
















AlO 


= 


3. 




P ( 1 ) 


+ 


64. 


* P ( 5 ) 
















A1 1 




4. 




P ( 1 ) 


+ 


Q 


. * P ( 2 ) 




1 A. 


if 


H( 3 ) 


4- 


P ( 4 ) 




A1 2 


= 


4 . 




P(1 ) 


+ 


22“^ 


. * P ( 2 ) 




256. 


if 


P { 6 ) 


+ 


P ( 6 ) 




AlP 


= 


4. 




P(1 ) 


+ 


122'=^ 


. * P(2) 


+ 


12P6. 


if 


P ( p ) 


4- 


P ( P ) 




A14 


= 


4 • 




P ( 1 ) 




25 


. * P ( 2 ) 


-f 


81. 


if 


P( 4 ) 


+ 


16. * 


P( 5 ) 


A1 5 


= 


4. 


if 


P( 1 ) 


+ 


441 


. PIP) 


+ 


62 A . 


if 


P( 4) 


+ 


16. * 


P ( 7 ) 


A16 


= 


4 • 


* 


P ( 1 ) 


+ 


4Q 


. P { 2 ) 


+ 


’56. 


if 


P( 5 ) 


+ 


81 . * 


P( 6 ) 


A17 




'S . 


if 


p n ) 




liiii 


. ^ P ( ’ ) 


+ 


226 . 


if 


P( 4) 


4- 


Q . * 


P(6 ) 




r «i| I 

iti 

m I 



, if <- >• ••rpiN ^ 

.^iji •• 

- I •- • •• ” 

^ • .tfi • 

■iw— I* • ' *0 * ** 



1 ^ . I •► • . 

l» » Al*- 
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DO B'^l 


1 = 1,4 


12 = 2 


*1 + 1 


Y2 = 2 


* I - 1 


Bl ( I ) = 


256. * P(l) * PN0RM/(Y2 * 3 , i 4 1 6*-*6 ) 


P2 ( I ) = 


4./(3. * (1. - U**2) * P( 12) ) 


301 B3 ( I ) : 


= -(16. * Y?**2 )/( 3. 1416**2 ) 


RETURN 




FND 




C 


SUBROUTINE FOR COMPUTING THE COEFFICIENTS 


C 


FOR THE USE OF THE CR TUT REDOX METHOD 


» LABEL 




* LISTS 





* FORTRAN 

SUBROUTINE SESAME 

COMMON W, AX, Bl, B2 , B3 ’ Al, A2 , A3, A4, A5, A^ , A7 , Ag, A<^, A10» 
1 All, A12, A13, AT 4, AIS, A] 6, A] 7, Mr,, A, C 

DIMENSION W(4), Bl(4), B2(4), B3(4), A(4,4), C(4) 

GO TO ( 401 , 400) 5 MG 



400 A ( 1 , 1 ) = 


: B 3 ( 1 ) *W ( 1 ) 


A(2 , 1 ) = 


:B3 ( 2 ) *W ( 2 ) 


A(3»l ) = 


:83 ( 3 ) *W ( 3 ) 


A (4, 1 ) = 


= B3 ( 4 ) *W (/..) 



GO TO 402 



4 O 1 A ( 1 , I ) = 


:B2(1)+B7(1)*AX+7,*A2*W(1)**2-a.*A1*W(1)*W(?)+A11*W(7)**7 


1 


+ A I 4*W ( 7 ) ** 2 +A 16 *W ( 4 )**7-7 .*A« + W ( 2 ' *W ( 7 ) -2 .*A10*V, ( 7 ) *N ( 4 ) 


A ( 2 , 1 ) = 


:-3.*Al*W(l )**2 + 2.*A11*va/(1 )*.V(2)-2,*AR*W(1)*W(3) 


1 


+2 A6*W ( 7 ) *W ( ^ ) -2 . *AO*W ( 2 ) ( 4 ) +A1 7*W ( 3 ) *W ( 4 ) 



A ( 7 , 1 ) = - 7 , *A q*W ( I ) *W ( ? ) + 2 . *AI 4*W ( I 1 *',v ( 7 ) - / . #A 10*W ( I ) u ) 
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+ A6*W(? )**2 + A17*W(2 ) *W ( A ) 



A(A«1)=-2.*A10*W(1)*W(^)+7,*a16*W(1)*W(A)-AQ*W(2)**2+A17*W(2)*W(^) 
402 A(1»2)=-^.*A1*W(1)**2 + 2.*A11'»^W(1)*W(?)+2.'^A6*W(2)*W(3) 

1 -2.*AO*W( 2 )*W( 4)-2.*A8*W( 1 )*W( 3 )+A17*W ( 2 )*W(4 ) 

A(2»2)=B2<2)+B^(2 )*AX + A1 1* W ( 1 ) ** 2 + ^ . * A3 *W ( 2 ) ** 2 + A 1 2* W ( 3 ) ** 2 
1 +A1 5*W ( 4 ) ^f*2 + 2 .*A6*W( 1 ) *W ( 3 ) ~2 . *AQ*W ( i )*W ( 4 ) 

A(^»2)=-AR*W(l )**2 + 2.*Aft*W( 1)*W(2)+2,*A12*W(2 )*i^'(3) 

1 +A 17*W( 1 ' *W ( 4 ) +2 . *A7*W ( 3 ) *W ( 4 ) 

A(4»2)=-2.*AQ*W(1)*W(2)+2.*Al3*W(2)*W(4)+A7^fW(2)**2 + Al7*>*Mi)«-W(-5) 

A ( 1 , 3 ) = A6*W ( 2 ) **2 + 2 .*A14*W ( 1 ) «-W( ^ ) -2, *A8*W ( 1 ) ^f-V/ ( 2 ) 

1 -2 . *A10*W ( ] ) *W ( 4 ) +A1 7*W ( 2 ) *W ( 4 ) 

A ( 2 » 7 ) = -AR*W ( 1 ) **2 + 2 . *Ai 2*W ( 2 ) *W ( ^ ) + 2 .*A7-- W ( 3 ) *W ( 4 ) 

1 +2.*A6*W(1)*W(2)+A17*W(1)*W(4) 

A(^»3)=B2(7)+B'3(-5)*AX + A14*W(1)**2+A12*W(2)**2 + 3,*A4*W(3)**2 
1 +A13*W(4)**2+2.*A7*W(2 )*W(4) 

A < 4 » 3 ) = -A 1 0* W ( 1 ) **2 + 2 .*A7*W ( 2 ) *W ( "s ) + 2 , *Al 3*W ( 2 ) *W ( 4 ) + A 1 7*W ( 1 ) *W ( 2 ) 
A(1 »4)=-A0*W(2)**2 + 2.*A16*W(1 )*W(4)-2.*A10*W( 1 )*W(3)+A17*W(2)*‘‘-'(3) 
A ( 2 > 4 ) = A7*W ( 3 ) **2 + 2 . *A1 R*ia/ ( 2 ) *W ( 4 )-2 .*AO*W ( 1 ) *W ( 2 ) + A1 7*W ( 1 ) *W ( 2 ) 
A(3,4)=-Ai0*W(1)**2+2.*A13*W(^)*W(4)+A17*W(1)*W(2)+2.*A7*W(2)*W(2) 

A(4»4)=B2(4)+R3(4)*AX+A]6*W(1)**2+A12*W(2)**2+A12#W(3)**2 
*1 +3.*A5*W( 4) **2 

c ( 1 ) = R1 (1 ) -R2 ( 1 )*W ( 1 ) -R2 ( 1 ) *w ( 1 ) *AX-A 2*W ( 1 ) **24-2 ,*A1*W ( 1 ) **2*W ( 2 ) 

1 -A n *W ( 1 ) *W( 2 ) **2-A6*W ( 2 ) #*2*W ( ^ ) -t-AO*W ( 2 ) **2*W ( 4 ) 

2 -A 1 4*W ( 1 ) *W ( ^ ) **2-A 16*W ( 1 ) *W ( 4 ) **2 + 2» *AR*W ( T ) *W ( 2 ) *, ( o ) 

3 +2 . *A 10*W ( 1 ) *W ( 3 ) *W ( 4 ) -A17*W ( 2 ) * W ( 3 ) *W ( 4 ) 

C ( 2 ) =B1 ( 2 ) -R 2 ( 2 )*W ( 2 )-R3 ( 2 ) *AX*W( 2 ) + A 1*W( 1 ) ** 2 - a 1 1 *W ( 1 ) **2*’v ( 2 ) 

1 -f A P*W ( ] ) **2*W ( 2 ) - A2 *W ( 2 ) * *2 - A ^ 2* ^^/ ( 2 ) * W ( 2 ) **2-A’’*W ( 2 ) **2*a ( 4 ) 

2 - A 1 3 * W ( 2 ) * ' ( 4 ) * * 2 - 2 . * A A * >' ( 1 ) * W ( 2 ) * W » 3 ) + 2 . * A c * W ( 1 )*.','( 2 ) * 4 ( 4 ) 
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00 



3 -A17*W ( 1 ) *W( 3 ) *W( 4 ) 

C(^)=R1(3)-B;?(3)*W('5) -B3 ( 3 ) *AX*W ( 3 ) + A ( 1 ) **2*»'‘ ( 7 ) 

1 -A14-«-W(l)**2*W(?)+AlO*W(l )«•*?•(■ W(A)-A6*W(1)*W(7)*^<-2 

2 -A 1 2«-W ( 2 ) **2*W ( 3 ) -a4*W ( R ) -K-JS-R-A 1 *W ( ^ ) *W ( 4 ) ■»^* 2 

R -A17«-W ( 1 ) *W( 2 ) *W ( 4 ) -2 .*A7^*-W ( 2 ) *W ( 3 ) *W ( A ) 

C(4)=P1 (4)-R2(4)*W(4)-R'»(4)^i•AX^(•W(4)+Al0*W<] )**2*W(3) 

1 -A16*W ( 1 ) ( 4 ) +A9*W ( 1 ) *W ( 2 ) **7-A ] 5*W( 2 ) ( 4 ) 

2 -A''*W(2)^^W(3)*-i^2-Al3*W(3)^*2*W(A)-A5*^W(A)**3 

3 -A 17*W ( 1 ) *W ( 2 ) ( 3 ) 

RETURN 

END 
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M3888-3401 »FMS»DFBUG*1*1 »500»0 

SOLUTION OF THE OEFLECTION EQUATION 

+ V sv-nlj- t V3^3 Sffl ^ sm -2M 

+ \J, 5 si-n ^ +Vl5^, si-n 

THIS METHOD IS S'^T UP 5=^0? A CONTINUOUS 
CHANGE IN AXIAL LOAD 'FROM ZERO TO A VALUE WITHIN 
THE VICINITY OF THE CRITICAL BUCKLING LOAD. ZERO 
LATERAL LOAD IS EXCLUDED. 

MAIN PROGRAM 

XEQ 

LABEL 

LIST8 

FORTRAN 

READ 1, IND 

1 FORMAT (15) 

PRINT 2 

2 FORMAT (114H1 PNORM AXIAL W(l9l)/H W(lt3)/H W(3»l)/ 

IH W(3»3)/H W(1.5)/H W(5.1)/H STRAIN ERROR ) 

DO 44 JOB = 1 » I ND 

READ 3» U. ALPHA, PNORM, LOOP, LAP 

3 FORMAT (3F8. 0,218) 

DIMENSION W(6), R(54), Bl(4), R2(5), B3(3), A(6,6)’ C(6), X(6) 

CALL CONST (U, ALPHA, PNORM, R, Bl, B2, B3 ) 

INITIAL ESTIMATES OF W ( I ) ' S 
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W( 1 ) 


= 0.100 


W( 2 ) 


= 0.010 


W( 3 ) 


= 0.010 


W ( 4 ) 


= 0.001 


W( 5 ) 


= 0.001 


W( 6 ) 


= 0.001 


C 


THE following IS THE ROUTINE FOR 


C 


CALCULATING THE DEFLECTION COEFFICIENTS FOR 


r 


chosen values of AXIAu load and NORMAL LOAD 



DIMENSION Z(3) » S(3»6) 



AX n 


0. 


MG = 


1 


RONE 


= R(1 ) 



DO 10 LL = 1,3 

4 CALL sesame 

CALL CROUT (A, C, X, M) 

GO TO ( 5 , 42 ) , M 

5 ERROR = 0. 

DO 6 J = 1 ,6 

6 ERROR = ERROR + ABSF(X(J)) 

IF (ERROR - 0.000001) 9, 9, 7 

7 DO 8 I = 1,6 

8 W ( I ) = W ( I ) + X ( I ) 

GO TO 4 



9 CALL 


CPRNT (PNORM, AX, W, ERROR, '^ON'^) 


CALL 


STORE (AX, W, Z, S, LL ) 


AX = 


AX + 0.5 



10 CONTINUE 
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DFL = 0.5 
KO = 1 

K< = 0 

11 DO 12 1=1,6 

FI = (5(2,1) - S(1,I))/(Z(2) - Z(D) 

F? = (5(3,1) - 5(2,1 ) )/(Z(3) - Z(2)) 

F3 = ( Z ( 3 ) - Z ( : ) ) /2. 

F4 = ( F2 - FI ) /F3 

12 W(I) = W(I) + F2 * DEL + 0.5 * F^ * DEL**2 
CALL EXCH (Z, 5) 

13 CALL 5F5AME 

CALL CROUT (A, C, X, MJ 
GO TO (14, 42 ) , M 

14 ERROR = 0. 

DO 1 5 J = 1 ,6 

15 ERROR = ERROR + AB5F(X(J)) 

IF (ERROR - 0.000001) 19, 10, 16 

16 KK = KK + 1 

IF (KK - LOOP) 17, 17, 40 

17 DO 18 1=1,6 

18 W( I ) = W( I ) + X( I ) 

GO TO 13 

19 IF ( W n ) - 5( 2 ,1 ) ) 22 » 22 , 20 

.20 CALL CPRNT (PNOR^S AX, W, ERROR, RONE) 

Z ( 3 ) = AX 

DO 21 1=1,6 

21 5(3,1) = W ( I ) 

AX = AX + DEL 
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KK = 


0 








GO TO 


1 1 






22 


GO TO 


( 23 . 


25 * 


27 ) 


23 


DEL = 


0.05 








DO 24 


I = 


1 » 6 






W( I ) 


= S(2 


♦ I ) 




24 


S(3» I 


) = W( I ) 






AX = 


Z( 2) 








Z(3) 


= AX 








AX = 


AX + 


DEL 






KO = 


2 








KK = 


0 








CALL 


EXCH 


( Z, 


S) 




GO TO 


1 13 






25 


DEL = 


0.01 








DO 26 


1 I = 


1 .6 






W( I ) 


= S(2 


♦ I ) 




26 


S(3, I 


) = W( I ) 






AX = 


Z ( 2 ) 








Z(3) 


= AX 








AX = 


AX + 


DEL 






KO = 


■x 








KK 


0 








CALL 


EXCH 


(Z» 


S) 




GO TO 13 







THE FOLLOWING IS ThE ROUTINE WHEN IN THE 
VICINITY OF THE BUCKLING LOAD. W(l) IS f-‘ADE 
The independent VARIABL'^ in place of the axial 
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LOAD 



27 DO 28 I = 1,6 
S(3,I) = 5(2,1) 
5(2,1) = 5(1,1) 

28 W( I ) = 5(3,1) 
Z(3) = Z (2 ) 

Z ( 2 ) = Z ( 1 ) 

AX = Z( 3 ) 



THP FOLLOWING VALUF5 O'" DLl AND DL2 ARE ARBITRARY. 
THE FIR5T VALUE OF DLl I5 SO CH05EN 5UCH THAT WHEN 
W(l) 15 DECREA5ED BY DI2* THc CON5EOUENT ITERATED 
50LUTI0N WILL NOT GO BACK ON THE PREVIOU5LY COYRUTED 
PATH. (NOTE - THE W(l) OR THE FIRST HARMONIC STARTS 
TO DECREASE IN VALUE BEFORE THE CRITICAL BUCKLING 
LOAD IS REACHED) 



DLl = 0.20 
DL2 = O.OlO 
MG = 2 



THIS LAST DO LOOP IS HOPED TO GO BEYOND THE CRI 
TICAL RANGE. START WITH LAP = 30. 

1 ,LAP 



DO 39 LAST 

KK = D 

AX = AX + DLl 
W( 1 ) = W( 1 ) - DL? 

DO 2'^ J = 2,6 

29 W(I) = W(I) + (S(3,J) - 5(2, J)) 
CALL EXCH (Z, S) 



30 CALL SESAMF 
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CALL CROUT (A, C, X, M) 



GO TO ( 31 , 42 ) » M 

31 «^RROR = 0. 

DO 32 J = 1,6 

32 ERROR = ERROR + ARSF(X(J)) 

IF (ERROR - 0.000001) 36, 36, 33 

33 <K = KK + 1 

IF (KK - LOOP) 34, 34, 40 

34 DO 35 I = 2,6 

35 W(I) = W(I) + X(I) 

AX = AX + X ( 1 ) 

GO TO 30 

36 CALL CPRNT (PNORM, AX, W, ERROR, R )NE ) 

IF (AX - Z( 2 ) ) 44, 44, 37 

37 Z(3) = AX 

DO 3R 1=1,6 

38 S( 3 , I ) = W( I ) 

DLl = Z (3 ) - Z (2 ) 

39 CONTINUE 
GO TO 44 

40 PRINT 41 

41 FORMAT (40H INFINITE LOOP OR NEEDS MORE ITERATION) 
CALL CPRNT (PNORM, AX, W, ERROR, RONE) 

GO TO 44 

42 PRINT 43 

43 FORMAT (19H CROUT WENT CRAZY) 

44 CONTINUE 



CALL FXIT 
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COMMON W, AX, Bl, B2» B3 » R» MG, A, C 
END 



SUBROUT I NES 



SUBROUTINE FOR SOLVING THE LINEARIZED SYSTEM 
OF equations by The crout redox method 



100 

101 



102 



103 



104 

105 



LABEL 

LISTS 

FORTRAN 

subroutine crout (A, c, X, M) 

DIMENSION A(6»6)» C(6)^ AA(6»5)> CC(6)» X(6) 
DO 1 on 1 = 1,6 

AA( I »1 ) = A( I ,1 ) 

DO 101 J = 2»6 
AA(1,J) = A(1,J)/A(1,1) 

DO 102 I = 2,6 
DO 102 J = 2,6 
AA( I , J) = 0. 

J = 2 
I I = J 

DO 105 1=11,6 

LIMl = J - 1 
DO 104 K = l,LIMi 

AA(I,J) = AA(I,J) + (AA(I,<) * AA(K,J)) 

AA ( I , J) = A ( I , J ) - AA ( I , J ) 

IF ( AA( J,J) ) 106, 116, 106 



106 IE (6 



J) no, 110, 107 
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107 I = J 

J = J + 1 
JJ = J 

DO 109 J = JJ»6 
LIM2 =1-1 
DO 108 K = 1»LIM2 

lOR AA(I»J) = AA(I,J) + (AA(I»K) * AA(K»J)) 
AA( ! »J) = A( I »J) - AA( I , J) 

109 AA(I»J) = AA< I »J)/AA( I » I) 

J = I + 1 

GO TO 103 

110 DO 111 I = 1»6 

X ( I ) = 0. 

111 CC( I ) = 0. 

CC ( 1 ) = C( 1 ) /AA( 1 »1 ) 

DO 113 I = 2»6 
LIM3 =1-1 
DO 112 K = 1,LIK3 

112 CC(I) = CC(I) + (AA(I»K) CC(K)) 

CC( I ) = C ( I ) - CC( I ) 

113 CC( I) = CC( I )/AA( T * 1) 

X(6) = CC(6) 

DO 115 I = 1 »5 
11=6-1 
LIMA = I I + 1 
DO 114 K = LIM4,6 

114 X(II) = X(ll) + (AA(II»K) * X(K)) 

115 X( I I ) = CC( I I ) - X( I I ) 
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M = 1 

GO TO lie 

116 PRINT 117 

117 format (28H SINGULARITY - NO SOLUTION) 

M = 2 

118 RETURN 
END 

C SUBROUTINE FOR STORAGE OF THREE SETS OF VALUES 

C OF AXIAL load AND THE FOUR DEFLECTION COEFFICIENTS. 

* LABEL 

* LIST8 

* FORTRAN 

SUBROUTINE STORE (AX, W, Z, S, LL ) 

DIMENSION W(6), Z(3), S(3,6) 

Z (LL ) = AX 
DO 203 I = 1,6 
203 S(LL, I ) = W( I ) 

RETURN 

END 

C THIS subroutine IS JUST A SHIFTING PROCEDURE 

C OF VALUES FROM ONE STORAGE LOCATION TO ANOTHER. 



* LABEL 

* LIST8 

» FORTRAN 

SUBROUTINE EXCH (Z, S) 
DIMENSION Z ( 3 ) , S ( 3 ,6 ) 
DO 204 I = 1,6 



= 1,2 



DO 204 J 









m 



j 

I 







k 
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JA = J 


+ 1 




204 


S( J, I ) 


= S( JA, 


I ) 




DO 205 


IJ = 1 


,2 




JB = IJ 


+ 1 




205 


Z ( I J) = 


Z ( JB ) 






DO 206 


JI = 1 


» 6 


206 


5(3, JI ) 


= 0. 






Z(3) = 


0. 






RETURN 








END 







c SUBROUTINE FOR PRINTING 

* LABEL 

* LISTS 

* FORTRAN 

SUBROUTINE CPRNT (PNORM, AX, W, ERROR, RONE) 

DIMENSION W(6) 

SUM=W ( 1)^^*2+W(2)**2+W(5)**2 + 9.*(W(3)-«-*? + W(4)**2)+25.*W(6)**2 
STRAIN = AX + 3.1A16**2 * SUM/(8. > RONE) 

PRINT 200, PNORM, AX, (W(I), 1=1,6), STRAIN, ERROR 

200 FORf'AT (F8.2, F11.5, 6F12.6, F11.5, F12.7) 

RETURN 

END 

C SUBROUTINE FOR COMPUTING ALL THE NECESSARY 

C CONSTANTS IN THE MAIN PROGRAM 

* LABEL 

* LISTS 

» FORTRAN 



SUBROUTINE CONST (U, ALPHA, PNORM, R, Rl, B2, B3) 
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DIMENSION P(13)» R(54), Bl(4)» B2 ( 5 ) , B30) 
P(l) = ALPHA**2 
P(2) = l./P(l) 

DO 300 I = 3»7 
Yl = (I - 2 ) **? 

300 P(I) = P(i)/(P(1) + Yl)**2 
DO 301 J = 8»11 

Y2 = (J - 6)**2 

301 P(J) = P(l)/(1. + Y2 * P(l))**2 
P(12) = P(1)/(A. * P(l) + 9.)**2 



P( 13 ) 


= 


P ( 1 ) / ( 4. + 


9 


, * P( 1 ) )**2 


R ( 1 ) 


= 


P ( 1 ) 






R ( 2 ) 


= 


P ( 2 ) 






R ( 3 ) 


= 


P ( 3 ) 






R ( 4 ) 


= 


P (4 ) 






R( 5 ) 


= 


P ( 8 ) 






R ( 6 ) 


= 


P ( 1 ) 




P ( 2) 


R ( 7 ) 


= 


P (1 ) 




81.^P(2) 


R(8 ) 


= 


P ( 1 ) 


+ 


625. *P ( 2 ) 


R(9) 


= 


81 .*P ( 1 ) 




P (2) 


R ( 10 ) 




625.*P(1 ) 


+ 


P ( 2 ) 


R ( 11 ) 


= 


P ( 1 ) 


+ 


P(8) 


R( 12 ) 


= 


P ( 1 ) 


+ 


81. *P (8) 


R( 13 ) 


= 


P (2 ) 




P ( 4) 


R( 14 ) 


= 


P (2 ) 




81.*P (4) 


R ( 18 ) 


r: 


P ( 1 ) 


+ 


3.*P (4) 


R ( 16 ) 


= 


P (2 ) 




3.*P ( 8) 


R( 17 ) 


= 


3.*P ( 1 ) 


+ 


64.*P (3) 
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f €iri** 


_ 46*1 

lim -• 


1 • 


U|* 

"V 


1 


•tr- 


• *0 



R ( 18 ) 


= 


3.*P (2 ) 


+ 


64.*P(3) 










R( 19) 


= 


64. *P ( •=> ) 




274.*P(4) 










R ( 20 ) 


r 


64. *P ( 3 ) 




274. *P ( 8 ) 










R( 21 ) 


= 


P (4) 


+ 


P ( 8 ) 










R(22 ) 


= 


P ( 3 ) 


+ 


D ( 5 ) 










R( 23 ) 


= 


P ( 3 ) 


-1- 


P ( 9 ) 










R ( 24 ) 


= 


9.*P ( 1 ) 




64 . *P ( 3 ) 


f 




25.*P ( 8 ) 




R (25 ) 


= 


9.<tP ( 2 ) 


+ 


64 .*P ( 3 ) 


+ 




25.*P(4) 




R(26) 


= 


25.*P ( 1 ) 


+ 


441 .*P ( 8 ) 


+ 




256.*P ( 9) 




R(27) 


= 


25.*P(2) 


+ 


441 .*p ( 4 ) 


K 




256. *P( 5) 




R(28 ) 


= 


272 .»P ( 3 ) 


+ 


625.*P(4) 


+ 




625. *-P ( 8 ) 




R( 29 ) 


= 


9.*P ( 1 ) 




64.*P( 3 ) 


+ 




1225. *P (4) 




R ( 30 ) 


= 


9.*P ( 2 ) 


+ 


64 .*P ( 3 ) 


+ 




1225. *P ( 8 ) 




R( 31 ) 


= 


25. *P ( 1 ) 


+ 


64.*P ( 3 ) 






9.*P( 10 ) 




R( 32 ) 


= 


25.*P(2) 


+ 


64. *P ( 3 ) 


-j. 




9.*P( 6) 




R( 33 ) 


= 


64.1S-P ( 3 ) 




441 .^fP ( 4 ) 






169. *P( 13 ) 




R ( 34 ) 


= 


64.*P ( 3 ) 


+ 


441 .*P ( 8 ) 


+ 




169. ^^P ( 12 ) 




R(35) 


= 


25.*P( 1 ) 




576. *P( 3) 


+ 




8281. ^^P( 13) 




R( 36 ) 


= 


25.*P(2) 


+ 


576. *P ( 3 ) 






8231. *P( 12 ) 




R(37 ) 


= 


1552. *P(3) 


+ 


28561 .*P ( 12 ) 


+ 


28561**P ( 13 ) 




R(38 ) 


= 


9.*P ( 1 ) 


+ 


4.*P(2) + 




16 


. *P ( 3 ) + 


P( 8 ) 


R ( 39 ) 


= 


4.*P ( 1 ) 


+ 


9.*P(2) + 




16 


.*P(3) + 


P ( 4 ) 


R ( 40 ) 


= 


36. *P ( 1 ) 




36.*P(2) + 


225 


,*P(^) + 225.*P(8) 


R(41 ) 


= 


25.*P(1 ) 


+ 


4 .*P ( 2 ) + 




81 


. *P( 8 ) + 


16.*P(P) 


R ( 42 ) 




4.^fP ( 1 ) 




25.*P(2) + 




81 


. «p ( A. ) + 


16.*P( 5 ) 


R ( 43 ) 


= 


36.*P ( 1 ) 


+ 


P ( 2 ) + 


144 


. *P ( 9 ) + 


9.*P ( 13 ) 


R ( 44 ) 


= 


P( 1 ) 


+ 


36 «*P ( 2 ) + 


144 


. *P ( 5 ) + 


9.<fP ( 12 ) 


R (45 ) 


= 


2.*P( 1 ) 


+ 


16.*P(3) + 




4Q 


.*P(m-) + 


9 .*P ( 8 ) 




I 



R ( 46 ) 




2 


.*P(2 ) 


+ 


16 .*P ( 3 ) 




9.*P(4 ) 


+ 


49.*P( 8 ) 


R(47) 


= 


16 


. *P ( 3 ) 


+ 


81 .*P (4 ) 


+ 


P(8) 


-f 


81.*P( 10) 


1 R(48) 


= 


16 


.*P(3) 


+ 


P(4) 




8l.*P(6) 


+ 


81 .*P( 8 ) 


R(49 ) 


= 


4 


.*P( 1 ) 




225 .*P ( 2 ) 




256.*P(3 ) 


- 1 - 


P( 6) 


R(50) 


_ = 


225 


.*P( 1 ) 


+ 


4.*P(2 ) 


+ 


256.*P ( 3 ) 


-f- 


P( 10 ) 


R(51 ) 


= 


256 


. *P ( 3 ) 


+ 


2401 .*P (4 ) 


+ 


4096. *P(5 ) 


+ 


2401 .*P( 12 ) 


R(52 ) 


= 


256 


. *P ( 3 ) 




2401 .*P ( 8 ) 


+ 


4096. *P(9) 




2401. *P( 13 ) 


R(53) 


= 


2. 


*P ( U 


+ 


64.*P(3) + 


256 


.*P(5) + 225 


.*P ( 8 ) + 49. 


R(54) 


= 


2. 


*P( 2 ) 


+ 


64.*P(3) + 


225 


.*P(4) + 256 


.*P( 9 ) + 49. 



DO 302 I = 1.3 
Y3 = 2 * I - 1 

BUI) = 256. * P(]) * PN0RM/(Y3 * 3 . 1 41 ) 

302 B3(I) = -(16. * Y3**2 )/( 3. 1416**2 ) 

R1 (4 ) = B1 ( 2 ) /3. 

DO 303 J = 1,5 
II = 2 * J + 1 

303 B2(J) = 4./(3. * (1. - U**2) * P(I1)) 

RETURN 

END 

C SUBROUTINE FOR COMPUTING THE COEFFICIENTS 

C FOR THE USE OF THE CROUT REDOX METHOD 

* LABEL 

* LISTS 

» FORTRAN 

SUBROUTINE SESAME 

COMMON W. AX, Bl, B2, B3, R, MG, A, C 

DIMENSION W(6), Bl(4), B2(5)» 63(3), R(54), A(6,6), C(6) 



*P (12 ) 
*P ( 13 ) 



GO TO ( 401 , 400 ) , mg 
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400 A ( 1 . 1 ) : 


=B3 ( 1 ) #W ( 1 ) 


A(2»l) = 


=B3 ( 1 ) *W (2 ) 


A(3.1) = 


=B3 ( 2 )*W( 3 ) 


A( 4, 1 ) : 


=B3 ( 2 ) *W ( 4) 


A(5»l ) = 


=B3( 1)*W(5) 


A(6. 1 ) = 


=83(3)*W(6) 



GO TO 402 



401 A( 1 » 1)^ 


= B2 ( 1 ) + B3 ( 1 ) *AX + 3. *R ( 6 ) *W( 1 ) ( 2 ) *^ ( 1 ) *W( 2 ) 


1 


-6.*R(1)*W(1)*W(3)+R(38)*W(2) **2+R ( 3 9 ) ( 3 ) **2 


2 


+ 32 .MR ( 3 )*W ( 2 ) ( 3 ) + 81 .*R ( 21 ) *W ( 4) 2 - 18 .*R ( 1 1 ) ( 2 ) *W ( 4 ) 


3 


-18. *R ( 13 )*W ( 3 ) *W ( 4 ) -6.*R ( 16 ) «-W ( 2 )*W( 5 )-6.*R ( 15 ) 3 )*W ( 6 ) 


4 


+1 62 . *R ( 5 )*W ( 4 ) *W ( 5 ) +1 62 .*R ( 4) *W ( 4) *W ( 6 ) +R ( 41 ) *W ( 5 ) **2 


5 


+R ( 42 ) *W ( 6 ) **2 


A ( 2 > 1 ) = 


= -3.*R ( 2 )-^W n ) **2 + 2. "+R ( 38 ) *W ( 1 ) *W ( 2 ) +32 .*R ( 3 )*W( 1 ) *W( 3 ) 


1 


-18 .*R ( 1 1 ) *W ( 1 ) *W( 4 ) -R ( 25 ) *W ( -5 ) **2-2 .-*R ( 24 ) *W ( 2 ) ( ^ ) 


? 


+R ( 4 O )*W( 3 )*W ( A ) -6. *R ( 16 ) *W ( 1 ) *W ( 5 ; +2. *R ( 18 ) *W( 2 ) *W ( 5 ) 


3 


+R(20)*W(3)*W(5)+R(19)*W(3)*W(6)-144.*R(23)*W(4)*W(5) 


4 


-9.*R(45)*W(4)*W(6)-R(27)*W(6)**2 


A ( 3 ♦ 1 ) ^ 


= -3.*R(1)*W(1)**2 + 32.'*R(3)*W(1)*W(2)+2.^R(39)*W(1)*W(3) 


1 


-18.*R(13)*W(1)*W(4)-R(24)*W(2)**2-2.*R(25)*w(2)*»v(3) 


2 


+R ( 40 )*W ( 2 ) *W( 4 ) -6.*R ( 15 ) *-W ( 1 ) *W(6) +R ( 20 )*W ( 2 ) *W ( 5 ) 


3 


+ R ( 19 )*W( 2 )*W ( 6 ) +2.*R ( 17 )*'^’ f 3 ) *W ( 6 ) -9.*R (46 ) *W( 4 ) *W( 5 ) 


4 


-1A4.*R(22)*W(4)*W(6)-R(26)*W(5)**2 



A(4»l)=-lf3.MR(ll)*W(l)*W(2)-18.*R(13)*W(l)*vV(3) 



1 


+ 1 62 . *R ( 21 ) *W ( 1 ) *W( 4 )+R (40 )*W( 2 ) ( 3 ) +16? .*R ( 5 ) *»'. ( 1 ) *^^ ( 5 ) 


2 


+162.*R ( 4 )*W( 1 ) *W ( 6 ) -144.*R( 23 ) *W( ^ ) *W ( 5 ) 


3 


-o.*R ( 45 )*W ( 2 ) *W ( 6) -9.*R ( 46 )*W ( 3 ) *W ( 5 ) -144.*R ( 22 ) *W ( 3 ) */v ( 6 ) 



4 



+ 18. *R ( 14) *W ( 4) *W( 5 ) +18.*R ( 12 ) -f-W ( 4) *W ( 6 ) 
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5 


+72 9 .*R ( ;:i ) *W( 5 ) *W( 6 ) 


A ( 5 ♦ 1 ) : 


=-6.*R(16)^W(1)*W(2)+16?.*R(3)*W(i)*W(4)+R(]8)*W(?)**2 


1 


+R(20)*W(2)^^W(3)+9.*R( 14 ) (4) **2-144. ^^R ( 2 3 ) * W ( 2 ) *W ( 4 ) 


2 


-9.*R(46)*W(3)*W(4)+2.*R(41)*W(1)*W(5)-2,*9(26)*V. (3)*W(5) 


3 


+729.*R(21 )*W(4)*W(6 ) 


A ( 6 » 1 ) = 


= -6.*R ( 1 5 )*W ( 1 ) *'•/( 3 ) +162.*R ( 4 )*W ( 1 )*W ( 4 ) +R( 17 ) *W( 3 ) **2 


1 


+R(19)*W(2)*W(3)+9.*R(12)*W(4) **2-v .*R ( 45 ) *W ( 2 ) *W ( 4 ) 


2 


-144 .*R ( 2? ) *W ( 3 ) *W ( 4 ) + 2 .*R ( 42 ) *W ( 1 ) *W ( 6 ) -2 . *R ( 27 ) *W ( 2 ) *-'^ ( 6 ) 


3 


+729.*R(21 )*W(4)*W( 5) 


402 A( 1 »2 ) = 


:-3 .*R ( 2 ) *W( 1 ) **2+2. *R ( 38 ) *W ( 1 ) *W ( 2 ) +32 .*R ( 3 )*W( 1 ) *W ( 3 ) 


1 


-18.*R(11)*W(1)*W(4)-2.*R(24)*W(2)*W(3)-R(25)*W(3)**2 


2 


+R ( 40 )*W( 3 ) *W( 4)-6. *R ( 16 ) *W( 1 ) *W ( 5 ) +2 .*R ( 18 ) *W( 2 )*W ( 5 ) 


3 


+R(20)*W(3)*W(5)+R(19)*W(3)*W(6)-144.*R(23)*W(4)*W(5) 


4 


-9.*R(45)*W(4)*W(6)-R(27)*W(6)**2 


A( 2 »2 ): 


= B2 ( 4 ) +83 ( 1 )*AX + R ( 38 ) *W ( 1 ) **2-2 .*R ( 24 ) *W ( 1 ) *W( 3 ) 


1 


+ 3.*R ( 9 )*W( 2 ) **2 + R ( 28 ) *W( 3 )**?-486.*R ( 1 ) *W( 2 ) *W( 4 ) 


2 


+ R ( 43 ) *W( 4 ) **2 + 2. *R ( 1 8 ) *W ( 1 ) *W ( 5 ) -51 2. *R ( ) *W( 3 ) *W ( 5 ) 


3 


-2.*R ( 29 )*W ( 3 ) *W( 6 ) +R ( 50 )*W ( 5 ) **2 + R( 51 )*W( 6 ) **2 


A(3»2) = 


=16. *R ( 3 ) *W ( 1 ) **2-2. *R ( 24 ) *W ( 1 ) *W ( 2 ) -2 . *R ( 25 )*W ( 1 ) *W ( 3 ) 


1 


+ R(40)*W'1)*W(4)+2.*R(28)*W(2)*W(-2)+R(20)*/J(1)*/h(5) 


2 


+R ( 19 )*W( 1 )*W ( 6 )-512 .*R ( 3 ) *W(2 ) *W( 5 ) -2 .*R ( 29 ) *W ( 2 ) *W ( 6 ) 


3 


-2 . *R ( 30 )*W ( 3 ) *W ( 5 ) - 51 2 .*R ( 3 ) *W ( 3 ) •'<■.'/ ( 6 )+9.*R ( 54 ) *W ( 4 ) *W ( 5 ) 


4 


+9.*R(53)*W(4)*W(6)+656.*R(3)*W(5)*^(6) 


A(4,2)^ 


= -9.*R ( 1 1 )*W ( 1 ) **2+R ( 40 ) *w ( 1 ) *W ( 3 ) -^43 . *R ( 1 ) *W ( 2 ) **2 


1 


+ 18.*R(43)*W(2)*W(4)-144.*R(2'^)*W(1)*w(5)-9.*R(31)*/.'(5)**2 


2 


-9.*R(45)*W(1)*Vv(6)+9,*R(54)*W(3)*W(5)+o.#R(53)*W(^)*W(6) 


3 


-9.*R ( 33 )*W ( 5 ) *W( 6 ) 


A ( 5» 2 ) ^ 


=-3. *R ( 1 6 ) *W ( 1 ) **2+2 .*R ( 18 ) * V ( 1 ) *.•< ( 2 ) +R ( 20 ) *.V( 1 ) *W ( 3 ) 
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1 



-144.*R(23)*W(1)*W(4)-512.*R(^)*W(2)*VH3)-P(3D)*-.’(3)**2 



2 


+ R ,*R ( 54 ) *W ( 3 ) *W ( 4 ) + 2 .*R ( 50 ) *W ( 2 ) ( 5 ) +6 56 .*R ( 3 ) *W ( 3 ) ( 6 ) 


3 


-18. *R ( 31 ) *W ( 4 ) ( 5 ) -Q.*R ( 33) *VJ( 4) •" W ( 6 )-R ( 36) ^W( 6 ) **2 


A ( 6 » 2 ) : 


= R(iq)*W(l)*W(3)-Q.*R(45)*W(|)-!(-W(4)-2,*R(79)'«-W(2)*''^<3) 


1 


-256.*R(3)*W(3)**2+9.*R(53)*W(3)^W<4)-2.*R(27)*a(1)*W(6) 


2 


+ 2.*R(51)*W(2)*W(6)+656.*R(3)*W(3)*W(5)-0.*R(33)«-W(4)<^W(5) 


3 


-2.*R ( 36 ) »W( 5 )*W ( 6 ) 


A( 1 »3 ) = 


:-3.*R(1)*W(1)**2+2.*R(39)*W(1)*W(3)+32.*R(3)*W(1)*W(2) 


1 


-] 8 . *R ( 13 ) *W ( 1 ) *W ( 4 ) -R ( 24 ) *W ( 2 ) ^^*2-2 .*R( 2 5 ) *W ( 2 ) *W ( 3 ) 


2 


+ R ( 4 O ) *W ( 2 ) *W ( 4 ) -6 . »R ( 1 5 ) *W ( 1 ) vfw ( 6 ' +2 . *R ( 1 7 ) ( 3 ) **-W( 6 ) 


3 


+ R(20)*W(?)*W(5)+R(]9)^^W(2)*i^^(6)-9.*R(46)*W(4)*W(5) 


A 


-144.*R ( 22 ) «-W ( 4} 6 )-R ( 26 ) '(•WC 5 ) **2 


A( 2»'^ ) = 


= 16.*R(3)*W(1)**2-2,*R(25)*v'M1)*/;(3i-2.*R(24)^^W(1)*W(2) 


1 


+ R ( 4 O ) -^W ( 1 ) *W ( 4 ) +2 , ^«•R ( 28 ) ( 2 ) *W ( 3 ) +R ( 20 ) ( 1 ) *W ( 5 ) 


2 


+ R(lQ)*W(l)*W(e)-2.*R(30)*'A(^)^^W(5)-512.*R(3)i'^W(^)">W(6) 


3 


-512.*R(3)-»-W(2)*W(5)-2.*R(29)i^W(2)*W(6)+9.*R(54)*W(A)*/;(5) 


4 


+ 9.^^R(53)*W(4)*W(6)+656.*R(3)*W(5)^<-'W(6) 


A(3>3) = 


= B2 ( 2 ) +83 ( 2 ) *AX + R ( 39 ) *W( 1 ) **2-2 .<<-R ( ■’5 ) 1 ) 2 ) 


1 


+R ( 28 )*W( 2 )**2+3. *R ( 7 ) *W( 3 )**2-486.*R ( 2 ) *W ( 3 ) *W( 4 ) 


2 


+q.*R(44)*W(4)**2+2.*R(17)*W(1)*W(6)-2.*R(30)*W(2)*W(5) 


3 


-5 1 2 . *R ( 3 ) *W ( 2 ) *W ( 6 ) +R ( 52 ) *W ( 5 ) ** 2 *R ( 49 ) *W ( 6 ) **2 


A(4»3) : 


= -9.*R ( 13 )*W ( 1 ) ** 2 +R ( 40 ) *'W ( 1 ) *W( 2 ) -24 3.*R ( 2 ) 3 ) **2 


1 


+ 18.MR ( 44 ) *W ( 3 ) *W(4 ) - 9 .*R ( 46) *W ( 1 )*W ( 5 )-9.*R ( 32 ) 6 )**2 


2 


-144.*R(22)*W(1)*W(6)+9.*R(54)*W(2)*W(5)+9.*R(53)*W(2)*W(5) 


3 


-9. *R ( 34 )*W ( 5 ) *W ( 6 ) 


A ( 5 » 3 ) : 


= R(20)*W(i)*W(2)-«5.*R(A6)*W(1)*W(4)-24^^.*R(3)*W(2)**2 


1 


-2 .*R ( ^0 )*W ( 2 ) *W ( 3 ) +Q .*R ( 54 ) *W ( 2 ) *'•' (4)-2.*R(’6)*«‘/(1)*W(5) 


2 


+ 2 . *R ( 52 ) *W ( ^ ) *W ( 5 ) + 6*^6 . *R ( ^ ) *Vn 7 ) --f /i ( 6 ) - ”. *R < 34 ) - ( A ) *W ( 6 ) 



c 



^-1 iMwaiM 

I iimiM ^1 —•»•■• -♦ 

i»»i ••••• |«VIM^*"** *•’“ 

• (n* ""um^ T*^ * * “■Jli^ * ^rWI^ * ^•••••••* 

I -I |||U.-C> i-**#^-' • •^SPPi'*"*^ • 

I ^ , «l) ••♦rt- •w^t l-**l •••■ ' 

«» '*’ ® "*** ® ® ♦ 'i*M 



IH»- 



.ni.». • — 



m jmm. -• 



L. 



• •- 1 


^iH 


t-gftllMIMv w’"*' 


•411 M 






• •■^IjiiKi IP 






■jmr^ 


«ltfft*lfl • * 




1 1 


M.I 


1 ari 1^*®^ 








I 



oir» X 



jt*i -' *• • * *•••*«•• 

» > •- il 4 i •••l»**-****<*^^*'*'^^®"****^^ ^ 



I I 




^CHf <i«*CWli I •• 



• « • H. t V I 

iii««ti np"» 



4 ■ P • 

■»- 



III 

• «ii«i i» I 



M 1 M*> '•* •—«•'»•• ■• '** ^ 

I I IM04» •»|.^l !*•> 



I • 

lll'i 

'll # ' I I '^^°*^rf*** * • ••Y.Hnn^i** 

f l•l»•|f‘^¥•• *••= 

\ I I 



• It 

i^lli» 4 ii 
• !• 'f■■^ 









-2 .*R ( 35 )*W ( 5 ) *W( 6 ) 



3 



A ( 6 » 3 ) : 


= -3 . *R ( 1 5 ) *W ( 1 ) **2+2 . *R ( 17 ) *'/i ( 1 ) ( 3 ) + R ( 19 ) ( 1 ) *W ( 2 ) 


1 


-lH4.*R(22)*W(l)*W(4)-R(29)*W(2)**2-512.*R(3)*rt(2)*W(3) 


2 


+ 9.*R(53)*W(2)*W(4)+2.*R(49)*W(3)*W(6)+6 56.*R(3)*.m(2)*W(5) 


3 


-9.*R (34 )*W( 4 ) *W ( 5 ) -18 .*R ( 32 )*W( ( 6 )-R ( 35 ) *vV ( 5 ) **2 


A ( 1 »4 ) ^ 


= 1 6 2 . MR ( 2 1 ) *W ( 1 ) *W ( 4 ) - 1 ft , *R ( 1 1 ) * W ( 1 ) * W ( 2 ) 


1 


-18 . *R ( 13 ) *W ( 1 ) *W ( 3 ) +R ( 40 ) *W ( 2 ) *'.V( 3 ) +1 62 .*R ( 5 ) *W ( 1 ) */J ( 5 ) 


2 


+ 1 6 ?.*R (4)*W( 1 ) *W(6 ) -144, *R ( 2ft )* Alt 2 ) *W ( 5 ) 


3 


-9,*R (45 )*W ( 2 ) *W ( 6 ) -9,*R ( 46 ) *W ( 3 )*(• ( 5 ) -144 ,*R ( 22 ) *W( 3 ) *W ( 6 ) 


4 


+ 18.*R(14)*W(4)*W(5)+18»*R(12)*>^(4)*W(6) 


5 


-729. *R ( 21 ) *W ( 5 ) *W( 6 ) 


A ( 2 »4 ) = 


=-9. *R ( 1 1 )*W ( 1 ) **2+R (40 ) *W ( 1 ) *W( 3 ) -243 ,*R ( 1 ) *W( 2 ) **2 


1 


+ 2. *R ( 43 ; *W ( 2 ) *W( 4) -14^‘ ,*R ( 23 ) *W( 1 ) * W( 5 ) -9 ,*R ( 45 ) *W( 1 )*W ( 6 ) 


2 


+9. *R ( 54 )*W ( 3 ) *W( 5 ) +9.*R ( 53 ) *W ( 3 ) *W( 6 ) -9,*R ( 31 )*W ( 5 ) **2 


3 


-9,*R ( 33 )*W ( 5 ) *W ( 6 ) 


A( 3 »4 ) ! 


= -9.*R ( 1 3 )*W ( 1 ) **2+R ( 40 ) *W ( 1 )*W ( 2 ) - ’43 , *R ( 2 ) *W ( 3 ) ^- <^2 


1 


+ 18 ,*R ( 44 )*W ( 3 ) *W ( 4 ) -9 ,*R ( 46 ) *W ( 1 ) *W ( 5 ) -9 ,*R ( 34) */J ( 5 ) *W ( 6 ) 


2 


-144,*R(22)*W(])*W(6)+Q.*R(34)*W(2)*W(5)+9.*R(53)*-'/(2)*‘^((6) 


3 


-9,*R(32 )*W(6 )**2 


A(4»4) : 


=81.*B2(l)+R3(2)*AX+fti.*R(21)*W(l)**2+9.*R(43)*W(2)**2 


1 


+ 243, *R (6 )*W( 4 ) **2 + 9 .*R (44 ) >W( 3 ) *^''2 + 18 ,*R ( 14) *W ( 1 ) *W( 5 ) 


2 


+ 1 8 ,MR ( 1 2 ) *W ( 1 ) *W ( 6 ) +8 1 ,*R ( 47) *W ( 5 ) **2 + 81 . *R ( 48 ) *W ( 6 ) **2 


-3 


+2592 ,*R ( 3 ) *W ( 5 ) *W( 6 ) 


A ( 5 »4 ) ^ 


= 8 1 ,*R ( 5 )*W( 1 ) **2+18 ,*R ( 14) *W( 1 ) *W( ^ ) +9 ,*R ( 54) *W( 2 ) *W( 3 ) 


1 


-144.*R(23)*W(l)*W(2)-9,*R(46)*W(l)*l^(3) 


2 


+729 , *R ( 2 1 ) *W ( 1 ) *W ( 6 )-16 ,*R ( 31 ) *W ( 2 ) *W ( 5 ) 


3 


-9.*R ( 33 ) *W ( 2 ) *W ( 6 ) -9. *R ( ft4 ) *W ( 3 ) *’•'/ ( 6 ) +162 ,*R ( 47 ) *W ( 4 ) *.\l ( 5 ) 



4 



+ 2 592 .*R ( 3 ) ( 4 ) *W ( 6 ) 
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A(6»4): 


:A1 , *R ( 4 ) *W( 1 ) **3+1 A ,*R ( 17 ) *W ( 1 ) *W ( 4 ) - Q .*R ( 4 3 ) * .V ( 1 )*•.'( 2 ) 


1 


-1 44 . *R ( 2 2 ) *W ( 1 ) *W ( 3 ) +9 . *R ( 53 ) *W ( 2 ) *W ( 3 ) 


2 


+729.*R(21)*W(1)*w(5)-9,*R(33)*W(2)»W(5)-9.*R(34)*W(3)*^(5) 


3 


-1 8 . *R ( 32 ) *W ( 3 ) *W ( 6 ) +2592 . *R ( ^ ) *W ( 4 ) *W ( 5 ) 


4 


+162.*R(48)*W(4)*W(6) 


Ad »3): 


= -6.*R(16)*W(1)*W(2)+162.*R(5)*W(1)*W(4)+2.*R(41)*.V(1)*/('(5) 


1 


+ R ( 18 )*W( 2 ) **2 + R ( 20 ) *W( 2 ) *Vv ( 3 ) -144.*R ( 23 ) *w'( 2 ) *W( 4 ) 


2 


-9.*R(46)*W(3)*W(4)+9.*R(14)*W(4)**2-2.*R(26)*W(3)*W(5) 


3 


+729. *R ( 21 )*W (4 ) *W( 6 ) 


A(2»5): 


= -3.*R(16)*W(1) **2+2 .*R ( 18 ) *W( 1 ) *V.'( 2 ) +R ( 20 ) *W( 1 )*W ( 3 ) 


1 


-144.*R(23)*W(1)*W(4)-R(30)*W(3) **^-512 .*R ( 3 ) *W ( 2 ) * W ( 3 ) 


2 


+ 9 . *R ( 54 )*W ( 3 ) *W ( 4) + 2 . *R ( 50 ) *W ( 2 ) *W ( 5 ) +6 56 . *R ( 3 ) *W ( 3 ) *VJ ( 6 ) 


3 


-18.*R(31)*W(4)*W(5)-9.*R(33)*W(4)*W(6)-R(36)*W(6)**2 


A ( 3 * 3 ) : 


:R(20)*W(1 )*W(2)-<=i.*R(46)*W(1)*W(4)-2.*R<26)*W( 1 )*W(5) 


1 


-2 56.*R ( 3 ) *W( 2 ) **2-2 .*R (30 )*W ( 2 ) *W( 3 ) + 9.*R ( 54) *>v ( 2 ) *W(4 ) 


2 


+2.*R(52)*W(3)*W(5)+656.*R(3)*W(2)*W(6)-9.*R(34)*W(4)*W(6) 


3 


-2 ,*R ( 35 )*W ( 5 ) *W ( 6 ) 


A ( 4 » 5 ) : 


= 81.*R(5)*W(1)**2-144.*R(?3)*W(])*W(7)-9,*R(46)*^(1)*'^((3) 


1 


+18. *R ( 14) *W ( 1 ) *W ( 4 ) +729.*R ( 21 ) *W (1 ) *W (6 ) 


2 


+ 9.*R(54)*W(2)*W(3)-18.*R(31)*W(2)*W(5)-9.*R(33)*/i((2)*W(6) 


3 


-9.*R (34 )-^W ( 3 ) *W( 6) +162.*R ( 47 ) *W (4) *W ( 5 ) 


4 


+2592.*R(3)*W(4)*W(6) 


A ( 5 » 5 ) : 


= B2(5)+B3(1)*AX + R(41)*W(1)**2-2."+R(26)*W(1)*W(3) 


1 


+ R(50)*W(2)**2 + R(52)*W(3)**2-18.*R(31)*'^<(2)*W(4) 


2 


+ 81 .*R ( 47 ) *W( 4 ) **2-2 .*R (35 ) *W( 3 ) *Wv 6 ) +3.*R ( 10 ) *W ( 5 ) **2 


3 


+R ( 37 ) *W ( 6 ) **2 


A ( 6 » 5 ) 


=729. *R ( 21 )*W ( 1 ) *W (4 ) +656.*R ( 3 ) *W ( 2 ) *W( 3 ) 


1 


-9 , *R ( 33 )*W ( 2 ) *W ( 4) -q, *R ( 34 ) *W ( 3 ) *W ( 4 ) +1 2'^6 .*R ( d *\: ( 4 ) **2 
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? -7.*R(^S)*W(^)*W(5)-:?»'*R(36)*W(7)«/.'(A)+?,«R(^7)*a(=7)*'^(5) 

A ( 1 » 6 ) =-6. *R ( 1 5 )*W ( 1 ) *W ( 3 ) +162 .*R ( 4 )*W ( 1 ) *W ( 4 ) + 2 .*R ( 4? ) 1 ) ( 6 ) 

1 +R ( 1 7 ) *W ( 3 ) **2 + R ( IR ) *W ( 2 ) *W ( 3 ) -9 .*R ( 46 ) VJ ( 2 ) +W ( 4 ) 

2 -144.*R(22)*W(3)*W(4)+9.*R(12)*W(4)**2-2.«R(27)*//(2)*/J(6) 

3 +729.*R(21 )*W(4)*W( 5 ) 

A(2»6)=R(l9)*W(i)*W(3)-g.*R(45)*W(i)*W(4)-2.*R(27)*W(i)*W(6) 

1 -2 56.*R(3)*W(3)**2-2.*R(29)*W(2)«-'Vi3)+9.«-R(5^)*W(3)*W(4) 

2 +2.*R(51)*W(2)*W(6)+656.*R(?)*W(-3)i^W(6)-9.*R(^'^)^W(4)*W(5) 

3 -2.^^R ( 36 )*W ( 5 )*W( 6 ) 

A(3»6)=-3.*R(15)*W(1)**2+R(19)*W(1)*W(2)+2.*R(17)^^W(1)<^W(3) 

1 -1 44.*R ( 22 ) *W ( ] ) *W( 4 )-R (29 )*W( 2 ) **2-51 2.*R ( 3 ) *W ( 2 ) *W ( 3 ) 

2 +9.*R(53)*W(2)*W(4)+2.*R(49)*W(3)*W(6)+656.*R(3)*W(2)*W(5) 

3 -18 . *R ( 32 ) *v\/ ( 4 ) *W ( 6 ) -9 .*R ( 3 '+) *W ( 4 ) *W ( 5 ) -R ( 35 ) *W ( 5 ) **2 
A(4»6)=ai.*R(4)*W(l)**2-9.*R(45)*W(l)*W(2)-144.*R(22)*W(l)*W(3) 

1 +18. *R ( 12 ) *W ( 1 ) *W ( 4 ) +729 .*R ( 21 ) *W ( 1 ) *K ( 5 ) 

2 +9.*R(53)*W(2)*W(3)-18.*R(32)*W(3)*W(6)-9.*R(33)*W(2)*W(5) 

3 -9.*R(34)*W(3)*W(5)+162.*R(48)*W(4)*W(6) 

4 +2592.*R(3)*W(4)*W(5) 

A(5»6 )=729.*R(21 )*W( 1 )*W(4)+6 56.*R( 3)*W(2)*W( 3) 

1 -9.*R(33)*W(2)*W(4)-9.*R(34)*W(3)*W(4)+1296.*R(3)*W(4)**2 

2 -2 .*R ( 36 )*W ( 2 ) *W ( 6 ) -2. *R ( 3 5 ) *W ( 3 ) *V~ ( 5 ) +2 .*R ( 37 ) *W ( 5 ) *W ( 6 ) 
A(6»6)=B2(3)+B3(3)*AX+R(42)*W(1)**2-2.*R(27)*W(1)*W(2) 

1 +R(51)*W(2)**2+R(49)*W(3)**2-18.*R'32)*W(3)*W(4) 

2 +81 . *R ( 48 ) *W ( 4 ) **2-2 .*R ( 36 )*W( 2 ) *W ( 5 ) +R ( 37 ) *W ( 5 ) **? 

3 +3. *R ( 8 ) *: W ( 6 ) **2 

01=B1(1 )-B2(l)*W(l)-B3(l )*AX*W(1)-R(6)*W(1)**3 

1 +3.*R ( 2 ) *W( 1 ) **2*W( 2 ) +3 .*R ( 1 )*W ( 1 ) **2*W( 3 ) -R ( 3 8 ) *W ( 1 ) * W ( 2 ) **2 

2 -R ( 39 ) *W ( 1 ) *W( 3 ) **2-32 . *R ( 3 ) *W ( 1 ) * \/ ( 2 ) *W ( 3 ) 
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-81 ( 2 1 ) *W ( 1 ) *W ( 4 ) »*2 + 18 . *R ( 1 1 ) *W ( 1 ) * W ( 2 )*W ( U ) 



4 -!-18.*R(13)*W(1)*W(3)*W(4)+R(24)*W(2)»^2*W(3) 

5 +R(25)*W(2)»W(3)**2-RU0)*W(2)*W(3)*W(A) 

6 +6.*R ( 16 )»W ( 1 ) *W( 2 ) *W ( 5 ) +6 . *R ( 15 )*.V ( 1 )*W ( ^ )*'a'( 6 ) 

7 -162.*R(5)*W(1)*W(4)*W(5)-162.*R(4)*W(1)*a/(4)*W(6) 

E1=-R(41)*W(1)*W(5)**2-R(4?)*W(1)*W(c,) * :^2-R ( 1 8 ) *W ( 2 ) ( 5 ) 

1 -R ( 17 ) *W ( 3 ) *^^2*W ( 6 ) -R ( 20 ) *W ( 2 ) *W ( ^ ) ( 6 ) 

2 -R ( 19 ) *W( 2 ) *W( 3 ) *W( 6 ) +144.*R ( 2 3 ) *W ( 2 ) ( 4 ) *W ( 5 ) 

3 +9 , *R ( 45 ) *W ( 2 ) ( 4 ) *W ( 6 ) +9 . *R ( 45 ) * W ( 3 ) *W ( 4 ) ( 5 ) 

4 +144 .*R ( 22 ) *W( 3 ) *W ( 4 ) *W ( 6 ) -9 .^^R ( 14 ) *W ( 4 ) **2*W ( 5 ) 

5 -9.*-R(12)*W(4)#*2*W(6)+R(27)*W(2)*W(6)**2 

6 +R ( 26 ) *W ( 3 ) *W ( 5 ) **2-729 .*R ( 2 1 ) *W ( 4 ) * W ( 3 ) * W ( 6 ) 
D2=Bl(2)-B2(4)*W(2)-B3(l)*AX*W(2)+R(2)*«'/(i)**3 

1 -R(38)*W(1)**2*W(2)-16.*R(3)*W(1)»*2*W(3) 

2 +9.*R(11D*W(1)**2*W(4)+R(25)*W(1)*W(3)-'>*2 

3 +2 .*R ( 24 ) *W ( I ) ^f■W ( 2 ) *W ( 3 ) -R ( 40 ) *W ( 1 ) #W ( 3 ) *W ( 4 ) 

4 -R ( 9 ) *W( 2 )**3-R ( 28 ) *W( 2 )*W ( 3 ) **2-R ( 43 )*W ( 2 )*W( 4 ) **2 

5 +243.*R(1)*W(2) **2*W ( 4 ) +3 . *R ( 1 6 ) *W ( 1 ) **2*W ( 5 ) 

6 -2 .*R ( 1 8D*W ( 1 ) *W ( 2 ) *W( 5 ) -R ( 20 )* \/( 1 ) *W ( 3 ) *W ( 5 ) 

7 -R ( 19 ) *W ( 1 ) *W ( 3 ) *W ( 6 ) +144. *R ( 23 ) *W ( 1 ) *W ( 4 ) *Vn 5 ) 

F2=9.*R(45)*W(1)*W(4)*W(6)+R(27)*W(1)*W(6)**2 

1 +R ( 30 ) *W ( 3 ) **2*W ( 5 ) +2 56 .*R ( 3 ) *W ( 3 ) **Z*W ( 6 ) 

2 +512.*R(3)*W(2)*W(3)*W(5)+2.*R(29)*W(2)*W(3)*W(6) 

3 -9 ,*R ( 54 ) *W ( 3 ) *W ( 4 ) *W ( 5 ) -Q . *R ( 53 ) *W ( 3 )*W ( 4 ) *W ( 6 ) 

4 -R ( 50 )*W ( 2 ) *W ( 5 ) **2-R ( 5 1 ) *W ( 2 ) *W ( 6 ) **2 

5 -6 56.*R(3)*W(3)*W(5)*'^(6)+9.*R(31)*W(4)*W(5)**2 

6 +9.*R(33)*W(4)*W(5)*W(6)+R(36)*W(5)*^(6)**2 



0^ =B1 ( 2 ) -R2 ( 2 ) *W( ) - R3 ( 2 )*AX^^W ( 3 ) +R ( 1 ) *W ( 1 ) -*;-*3 
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1 - 16 . *R ( 3 ) *w ( 1 ) ( ? ) -R ( 3Q ) *w ( 1 ) * >7*W f 3 ) 

2 +9.*R( 13 ) *W ( 1 ) **2*W( 4 )+R ( 24 ) *W ( 1 ) *W ( 7 )**2 

3 +2 .*R( 25 )*W ( 1 )*W ( 2 )-'‘W( 3 ) -R ( 40 ) * V( 1 ) *.v (2 ) *W ( 4 ) 

^ -R ( 28 )*W ( 2 ) **2«-W ( 3 ) -R ( 7 ) *W ( 3 ) **3+2 43 .*K ( 2 ) *^ ( 3 ) ""-*2*A' ( 4 ) 

5 -9,*R(44)*W(3)*W(4)**2+3.*R(15)*W( 1)**2*W(6) 

5 -R ( 20 )*W ( 1 ) *W ( 2 ) *W( 5 ) -R ( 19 ) *W( 1 )*W( 2 ) *W( 6 ) 

7 -2 ,*R ( 1 7 ) *W ( 1 ) *W ( 3 ) *W ( 6 ) +9 . *R ( 46 ) *W ( 1 ) ( 4 ) *W ( 5 ) 

E3 = ’44,*R(22)*W(1)*W(4)*W(6)+R(26)*W( 1 )^fW(5)**2 

1 +256.*R(3)*W(2)**2*W(5)+R(29)*W(2) *^^2*'*' ( 6 ) 

2 +2 ,<^R ( 30 ) *W (2 )*W ( 3 ) *W( 5 ) +5 12 .*R ( 3 ) *W ( 2 ) *W( 3 )*W ( 6 ) 

3 -q,*R(54)-k-W(2)*W(4)*W(5)-9.*R(53)*W(2)*W(4)*W(6) 

4 -R ( 52 ) *W ( 3 ) *W ( 5 ) **2-R (40 )*W ( 3 ) *W ( 6 ) **2 

5 -656.*R(3)*W(2)*W(5)*^(6)+9.*R(32)*W(4)*W(6)**2 

6 +9.*R ( 34 ) *W (4 )*W( 5 ) *W ( 6 ) +R ( 35 ) *W ( 5 ) **2*W( 6 ) 

D4 = Bl(4)-81.*B2(l)*W(4)-B3(2)*AXHi-W(4)+9,-)<-R(ll)*W(l)*«-2*W(2) 

1 +9,*R ( 13 ) *W ( 1 )^t*2*W ( 3 )-81 .*R ( 21 ) *W ( 1 ) **2*W ( 4 ) 

2 -R(40)*W(1)*'''(2)*W(3)+81.*R(1)*W(2)**3 + 81.*R(2)*W(3)**3 

3 -q,*R(43)*W(2)**2*W(4)-81.*R(6)*W(4)**3-o,*R(44)*W(3)**2*'''((4) 

4 -81 ,*R ( 5 ) *W ( 1 ) **2*W ( 5 ) -8 1 . *R ( 4 )*W( 1 ) **^*W ( 6 ) 

5 +144 ,*R ( 23 ) *W ( 1 ) + W( 2 ) *W ( 5 ) +o.*R ( 4=. ) *W ( 1 ) *W ( 2 )*l'J ( 6 ) 

6 +9. *R ( 46 ) *W ( 1 ) *W ( 3 ) *W ( 5 ) +144.*R ( 22 ) *W ( i ) *W ( 3 ) *W ( 6 ) 

F4 = -18. *R ( 14 ) *W ( ] ) + W ( 4 ) *W ( 5 ) -18 .*R ( 12 ) *W ( 1 ) *V\' ( 4 ) *W ( 6 ) 

1 -72Q .+R ( 2 1 ) *W ( 1 ) *W ( 5 ) *W ( 6 ) -9 , *R ( 54 ) +W ( 2 ) *W ( 3 ) *W ( 5 ) 

2 -Q,*R ( 53 ) *W (2 ! *W ( 3 ) *W(6 ) +9.*R ( 31 )+W( 2 )*W ( 5 ) **2 

3 +Q,^^R ( 32 ) *W ( 3 ) + W ( 6 ) **2 + 9 .*R ( 33 ) HV ( 2 ) * .V ( 5 ) * W ( 6 ) 

4 +9.*R ( 34 ) *W ( 3 )*W( 5 ) *W ( 6 ) -81 .*R ( 47 ) *W ( 4 ) *W( 5 ) **2 

5 -8 1 . *R ( 48 ) *W ( 4 ) *W ( 6 ) **2-2592 . *R ( 3 ) *W ( 4 ) *W ( 5 ) *W ( 6 ) 
D5=Bl(3)-B2(5)*W(5)-B3(l)*AX*Vi(5)+3.*R(16)*W(l)**2*»^(2) 
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I 

I _fl].*R( 5 )*W(l)** 2 *W( 4 )-R( 18 )^^W(l)*W( 2 )** 2 -R( 20 )*w(l)*ft( 2 )^^j( 3 ) 
I 2 - 9 .*R ( 14 ) *W ( 1 ) ^<-W ( 4 ) **2 + 144 . *R ( 23 ) ( 1 ) *W ( 2 ) ( 4 ) 

! 3 +9 ,*R ( 45 ) *W ( 1 ) *W ( 3 ) *W ( 4 ) +2 56 • *R ! 3 ) *W ( 2 ) ** 2 * 1 ^ ( 3 ) 

{ 4 +R ( 30 ) *W ( 2 ) *' ' ( 3 ) ** 2-9 .*R ( 5^0 *W ( 2 ) *W( 3 ) *W ( 4 ) 

5 -R ( 41 ) *W ( 1 ) **2*W ( 5 ) +2 .*R ( 2 6 ) *W ( 1 ) *W ( 3 ) *W ( 5 ) 

5 -729«*R(21) *W ( 1 ) *W ( 4 ) *W ( 6 ) -R ( 5 0 ) *W ( 2 ) ** 2 *W ( 5 ) 

E5 = -R ( 52 ) *W ( 3 ) **2*W ( 5 ) -6 56 . *R ( 3 ) * W ( 2 ) ( 3 » W ( 6 ) 

1 +i 8 ,*R( 31 )*W( 2 )*W( 4 )*W( 5 )+ 9 .*R( 33 )*W( 2 )*W( 4 )*W( 6 ) 

2 +9 , *R ( 34 ) *W ( 3 ) *W ( 4 ) * W ( 6 ) -8 1 .*R ( +7 ) *W ( 4 ) **2*W ( 5 ) 

3 - 1296 .*R( 3 )*W( 4 )** 2 *W( 6 )+R( 36 )*W( 2 )*W(r >)**2 

4 + 2 ,*R( 35 )*W( 3 )*W( 5 )*W( 6 )-R( 10 )*W( 5 )** 3 -R( 37 )*W( 5 )*W( 6)**2 
D 6 =B 1 ( 3 ) -B 2 ( 3 ) *W ( 6 ) -B 3 ( 3 ) *AX*W ( 6 ) + 3 .*R ( 1 5 ) ( 1 ) ** 2 *W ( 3 ) 

1 _ 81 .*R(i)*W(l)** 2 *W( 4 )-R( 17 )*W(l)*W( 3 )** 2 -R(l 9 )*W(l)*W( 2 )*W( 3 ) 

2 - 9 ,*R( 12 )*W( 1 )*W( 4 )** 2 + 9 .*R( 45 )*W( 1 )*W( 2 )*W( 4 ) 

3 + 1 44 , -if-R ( 2 2 ) *W ( 1 ) *W ( 3 ) *W ( 4 ) +R ( 29 ) *W ( 2 ) **2*W ( 3 ) 

4 + 256 .*R ( 3 )*W ( 2 ) *W ( 3 ) ** 2 - 9 . *R ( 5 3 ) *W ( 2 ) *W ( 3 ) *W ( 4 ) 

5 -R ( 42 ) *W ( 1 ) ** 2 *W ( 6 ) +2 «*R ( 27 ) *W ( 1 ) *W( 2 ) *W ( 6 ) 

6 -7 29 . *R ( 2 1 ) *W ( 1 ) *W ( 4 ) *W ( 5 ) -R ( 5 1 ) *W ( 2 ) **2*W ( 6 ) 

E 6 = -R ( 49 ) *W ( 3 ) ** 2 *W ( 6 ) - 656 .*R ( 3 )*W( 2 ) *'W ( 3 ) *W( 5 ) 

1 + 9 ,*R( 33 )*W( 2 )*W( 4 )*W( 5 )+ 9 .*R( 34 )*-W( 3 )*W( 4 )*W( 5 ) 

2 + 18 .*R( 32 )*W( 3 )*W( 4 )*W( 6 )- 1296 .*R( 3 )*W( 4 )** 2 *W( 5 ) 

3 -81 .*R ( 48 )*W ( 4) **2*W ( 6 ) +R ( 3 5 ) *W ( 3 ) *W ( 5 ) **2 

4 +2 .*R ( 36 ) *W ( 2 ) *W ( 5 ) *W ( 6 ) -R ( 37 ) *W ( 5 ) **2^W ( 6 ) -R ( P ) *'w ( 6 ) **^ 

C { 1 ) = D 1 + El 
C( 2 ) = D 2 + E 2 
C ( 3 ) = D 3 + E 3 
C ( 4 ) = D 4 + E 4 



C( 5) 



D5 + E5 



•f • 



IF-* 






g» LT t l*^*-***!* 

rnltm, l*V^ 

■•••••-. .*rl •••► 



^ ^ ^ ** lfc«*ffr***»*^ 

I ^tiMi*****^*^'^* * ••• *^ flO^^**** ** 

-^1 •••'» •«•••“-• <j'*'**» 

, rt***'* *"*** 

»*»I-»*U«**»**'*""*»*'-'*^" 

• "HfH «•»!*►*• '• "I 

• I Pi*** • 

*'■“**•" * •'** 

l..> «ti.l 

I I 



■ 

i 

* 

I# 

I 

I 

i 



% 






114 



C(6) = D6 + E6 

RETURN 

END 

DATA 
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